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EXECUTIVE  SUMMARY 


Details  of  the  barotropic  hydrodynamic  model,  ADCIRC-2DDI,  including  formulation  of  the 
governing  equations  and  their  numerical  solution  are  presented.  Validation  of  ADCIRC-2DDI  pre¬ 
dicted  tidal  heights  and  currents  in  the  North  Sea/English  Channel  is  undertaken  in  the  context  of 
the  Tidal  Flow  Forum  as  outlined  by  Werner  and  Lynch  (1988).  Sensitivity  studies  examine  the 
effects  of  nonlinear  mechanisms,  forms  of  the  tidal  forcing,  and  the  behavior  of  various  model 
parameters  on  the  computed  tidal  response.  Mesh  resolution  issues  are  also  investigated  in  the 
context  of  tidal  dynamics.  Throughout  this  report,  model  experiments  and  the  discussion  of  results 
are  focused  on  improving  understanding  of  the  ADCIRC-2DDI  hydrodynamic  model  and  its 
application  to  coastal  tide  prediction. 


COASTAL  TIDE  PREDICTION  USING  THE  ADCIRC-2DDI  HYDRODYNAMIC 
FINITE  ELEMENT  MODEL:  MODEL  VALIDATION  AND  SENSITIVITY 
ANALYSES  IN  THE  SOUTHERN  NORTH  SEA/ENGUSH  CHANNEL 


1.0  INTRODUCTION 

Numerical  modeling  has  become  an  essential  tool  for  assessing  the  hydrodynamics  of  continental 
margin  waters.  It  is  important  to  recognize  that  the  computed  response  of  these  waters  is  controlled 
by  the  various  components  that  make  up  the  model,  including  the  governing  equations,  boundary 
conditions,  forcing  functions,  discrete  equations,  grid  structure,  and  the  computational  domain 
itself.  The  more  that  is  understood  about  a  numerical  model,  how  its  various  components  influence 
computations  and  is  put  into  practice,  the  more  successful  the  model  will  be  in  representing 
hydrodynamic  processes  within  shallow  waters.  The  work  described  herein  focuses  on  the  numerical 
modeling  of  tidal  circulation  in  coastal  waters.  In  particular,  sensitivities  of  the  tidal  elevations  and 
currents  computed  by  the  ADCIRC-2DDI  hydrodynamic  model  simulator  are  investigated.  Important 
considerations  in  the  numerical  modeling  of  tidal  dynamics  are  the  capture  of  nonlinear  interac¬ 
tions,  the  specification  of  boundary  conditions,  and  the  adequacy  of  the  grid  resolution.  The  North 
Sea  Benchmark  (Werner  and  Lynch  1988)  provides  the  backdrop  for  this  sensitivity  study  and 
validation  of  tidal  predictions.  Results  included  herein  provide  guidelines  for  future  applications  of 
the  ADCIRC-2DDI  hydrodynamic  model  to  the  simulation  of  tidal  circulation  in  the  coastal  ocean. 

Tides  are  long  period  waves  generated  by  gravitational  forces  of  the  sun  and  the  moon  on  the 
ocean  waters  (Hendershott  1981).  The  period,  wavelength,  and  amplitude  characteristics  of  the  tide 
depend  on  geometric  properties  of  a  specific  water  body,  i.e.,  the  coastal  outline  and  bathymetric 
profile.  Tidally  generated  water  heights  and  currents  can  be  the  dominant  feature  in  a  coastal  area 
or  they  may  serve  as  background  circulation  that  contributes  in  an  influential  manner  to  overall 
coastal  dynamics.  Tidal  predictions  are  a  necessary  component  of  any  description  of  the  coastal 
environment  for  navigational  purposes,  coastal  fisheries,  and  military  operations,  to  name  a  few 
applications.  The  success  of  each  of  these  endeavors  rests  upon  accurate  prediction  of  the  tidal 
response  of  the  coastal  ocean  and,  consequently,  on  the  formulation  of  an  accurate  numerical 
hydrodynamic  model. 

Justification  for  selection  of  the  ADCIRC-2DDI  finite  element  (FE)  numerical  model  is  based 
on  accuracy,  mesh  flexibility,  and  computational  efficiency.  ADCIRC-2DDI  (Luettich  et  al.  1992) 
implements  the  generalized  wave-continuity  equation  (GWCE)  and  momentum  balance  equations 
for  which  accuracy  is  well  documented  with  respect  to  the  solution  of  various  shallow-water  problems 
(Foreman  1988;  Lynch  et  al.  1988;  Walters  1988;  Werner  and  Lynch  1989;  Walters  and  Werner 
1989;  Gray  1989;  Lynch  and  Werner  1991;  Luettich  et  al.  1992).  In  addition,  the  FE  formulation 
of  ADCIRC-2DDI  leads  to  tremendous  grid  flexibility,  allowing  easy  incorporation  of  coastline 
detail  and  nodal  densities  that  range  over  3  to  4  orders  of  magnitude.  The  wide  variation  in  nodal 
density  permits  significant  resolution  of  shoreline  geometry,  as  well  as  high  levels  of  refinement 
near  shallow,  coastal  areas  and  in  regions  of  rapid  bathymetric  change.  Moreover,  the  discrete 
problem  remains  well  within  computational  limits  despite  the  large  variation  in  nodal  density. 
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Following  a  detailed  presentation  of  the  numerical  hydrodynamic  model,  its  development,  equations, 
and  numerical  solution,  validation  of  predicted  tidal  heights  and  currents  in  the  North  Sea/English 
Channel  is  undertaken  within  the  context  of  the  Tidal  Flow  Forum  as  outlined  by  Werner  and  Lynch 
(1988).  Discussion  of  sensitivity  studies  examines  the  effects  and  behavior  of  nonlinear  terms  in  the 
mode]  equations,  forms  of  the  tidal  forcing,  and  a  variety  of  model  parameters.  Issues  relating  to 
mesh  resolution  are  investigated  within  the  context  of  tidal  dynamics.  A  final  statement  summarizes 
major  conclusions  and  observations  that  augment  our  understanding  of  tidal  dynamics  and  increase 
our  capabilities  in  applying  the  ADCIRC-2DDI  hydrodynamic  model  to  coastal  tide  prediction. 


2.0  THE  HYDRODYNAMIC  MODEL 

The  purpose  of  a  barotropic  hydrodynamic  model  is  to  characterize  the  important  flow  features 
of  surface  waters  driven  by  tides,  wind,  and  atmospheric  pressure  gradients.  The  complexity  of 
circulation  patterns  in  coastal  and  continental  shelf  margin  waters  necessitates  use  of  a  numerical 
hydrodynamic  model.  Such  models  compute  spatial  and  temporal  distributions  of  velocity  and  sea 
surface  elevation  from  which  circulation  patterns  can  be  inferred.  Success  of  the  computed  tidal 
response  in  resembling  the  actual  coastal  ocean  depends  in  part  on  selection  of  an  appropriate 
hydrodynamic  model.  An  optimal  model  is  one  that  captures  tidal-induced  dynamics,  yet  remains 
computationally  feasible. 


2.1  Criteria  for  Selection 

The  theoretical  basis  of  a  hydrodynamic  model  is  found  in  the  principles  of  mass  and  momentum 
conservation.  For  vertically  well-mixed  surface  waters  experiencing  tidal  and  atmospheric  forcing, 
the  flow  physics  are  described  by  the  shallow-water  equations,  a  depth-integrated  form  of  the 
conservation  laws  (Le  Mehaute  1976).  The  shallow-water  equations  have  been  used  successfully  by 
engineers  and  researchers  for  many  years  to  predict  tidal  circulation  (e.g.,  Reid  and  Whitaker  1981; 
Gray  and  Kinnmark  1983;  LeProvost  and  Vincent  1986;  Flather  1988;  Foreman  1988;  Baptista 
et  al.  1989;  Al-Rabeh  et  al.  1990;  Westerink  et  al.  1989,  1992a,  1994a,  1995;  Luettich  and  Westerink 
1995). 

Naturally,  formulation  of  the  hydrodynamic  model  and  the  solution  strategy  for  the  discrete 
problem  must  be  accurate.  Wavelength  and  phase  propagation  characteristics,  mass  conservation 
properties,  and  performance  in  idealized  test  cases  and  field  applications  are  all  used  to  evaluate 
model  accuracy. 

A  second  consideration  in  selection  of  a  model  is  efficiency  of  the  numerical  solution  algorithm. 
Meaningful  boundary  forcing  can  necessitate  the  use  of  large  model  domains  that  cover  the  continental 
shelf,  coastal  regions,  and  include  portions  of  the  deep  ocean  (Blain  et  al.  1994).  Furthermore,  the 
time  frame  for  typical  simulations  ranges  from  hours  to  months.  The  feasibility  of  a  model  that 
utilizes  large  domains  and  extended  simulation  periods  depends  upon  the  efficiency  of  the  numerical 
solution  algorithm.  Efficiency  is  achieved  by  minimizing  the  number  of  degrees  of  freedom  and  the 
number  of  operations  per  degree  of  freedom  at  each  computational  timestep. 

Lastly,  both  efficiency  and  accuracy  of  the  hydrodynamic  model  formulation  are  intertwined 
with  the  characteristics  of  the  discrete  mesh  and  its  relation  to  the  numerical  solution  algorithm. 
Grid  flexibility  allows  minimization  of  the  number  of  degrees  of  freedom,  while  simultaneously 
providing  the  localized  resolution  needed  for  accurate  model  predictions. 
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The  selection  criteria  for  any  hydrodynamic  model  is  summarized  as  follows:  governing  equations 
based  on  conservation  principles,  an  efficient  and  accurate  numerical  solution  algorithm,  and  a 
significant  level  of  grid  flexibility. 


2.2  Current  Numerical  Modeling  Approaches 

The  search  for  accurate  and  efficient  solution  algorithms  for  the  shallow-water  equations  has 
led  to  a  variety  of  equation  formulations  and  the  use  of  several  numerical  solution  techniques.  Two 
prevalent  numerical  discretization  strategies  utilize  the  finite  difference  and  FE  approaches.  Finite 
difference  methods  discretize  derivative  operators  contained  in  the  model  equations  using 
point  difference  expressions  (Lapidus  and  Pinder  1982),  while  FE  techniques  approximate  the 
solution  to  the  model  equations  through  the  use  of  interpolation  functions  (Celia  and  Gray  1992). 
A  review  of  several  FE  solution  techniques  for  the  shallow-water  equations  is  given  by  Lee  and 
Froehlich  (1986).  Westerink  and  Gray  (1991),  in  their  recent  review,  note  increasing  similarities 
between  finite  difference-  and  FE-based  models,  especially  with  regard  to  the  accurate  propagation 
of  short  wavelengths,  grid  flexibility,  and  complex  boundary  representation.  Further  comparison 
and  discussion  of  these  methods  is  given  by  Blain  and  McManus  (1998). 

Early  finite  difference  solutions  of  the  shallow-water  equations  were  quite  accurate,  owing  in 
part  to  implementation  of  a  staggered  grid  approach,  which  successfully  avoided  the  introduction 
of  artificial  short  waves  (Westerink  and  Gray  1991).  In  contrast,  early  FE  schemes  were  plagued 
by  spurious  oscillations  or  numerical  noise  due  to  the  introduction  of  2  •  Ax  wavelengths,  a  conse¬ 
quence  of  folded  dispersion  characteristics  (Lynch  1983;  Kolar  et  al.  1994b).  Note  that  2- Ax 
wavelengths  are  the  smallest  wavelengths  captured  by  a  mesh  having  a  minimum  spacing  Ax.  Early 
FE  algorithms  included  excessive  damping  mechanisms  to  counter  the  generation  of  spurious  modes 
and,  as  a  result,  yielded  inaccurate  solutions  (Gray  1982;  Gray  and  Kinnmark  1983).  It  wasn’t  until 
Lynch  and  Gray  (1979)  introduced  the  wave-continuity  equation  (WCE)  formulation  for  the 
shallow-water  equations  that  the  viability  of  FE  approaches  improved. 

The  WCE  formulation  simply  involves  a  rearrangement  of  the  shallow-water  equations  prior  to 
spatial  discretization.  Solutions  using  the  WCE  successfully  suppress  short  wavelengths  without 
resorting  to  nonphysical  dissipation.  Understanding  the  origin  and  behavior  of  spurious  oscillations 
within  FE  formulations  (Walters  and  Carey  1983),  together  with  extensive  numerical  testing  of  the 
WCE  formulation  (Kinnmark  and  Gray  1984,  1985;  Luettich  et  al.  1992;  Westerink  et  al.  1992b; 
Kolar  et  al.  1994  a,b)  has  resulted  in  very  accurate  WCE-based  FE  solutions  (Lynch  and  Gray  1979; 
Foreman  1983). 

In  addition  to  accuracy,  efficiency  of  a  numerical  algorithm  must  also  be  considered.  Representation 
of  complex  circulation  patterns  by  a  hydrodynamic  model  often  requires  highly  refined  grids  having 
very  small  nodal  spacings.  As  a  consequence,  the  stability  criteria  for  various  time-stepping  schemes 
dictate  use  of  a  very  small  time  interval  for  computations.  The  practicality  of  using  small  timesteps 
for  long  period  simulations  depends  on  the  efficiency  of  the  numerical  solution  algorithm.  Implicit 
time-stepping  schemes  are  generally  more  stable  and  allow  larger  timesteps,  but  result  in  time- 
dependent  matrices  that  require  reassembly  and  resolution  at  every  timestep.  This  procedure  is 
computationally  intensive.  Many  finite  difference  models  overcome  this  problem  by  implementing 
the  alternating  direction  implicit  (ADI)  solution  algorithm  that  reduces  a  two-dimensional  (2D) 
problem  to  a  sequence  of  one-dimensional  (ID)  problems,  thus  significantly  reducing  the  compu¬ 
tational  effort  (Leendertse  1987).  The  ADI  approach,  however,  cannot  be  applied  to  unstructured 
grid  algorithms  such  as  FE  solution  strategies. 


In  contrast  to  an  application  of  FEs  to  the  primitive  shallow-water  equation,  a  WCE-based  FE 
formulation  leads  to  decoupled  elevation  and  velocity  solutions  and  sparse,  symmetric  matrices 
(Lynch  and  Gray  1979).  Solutions  for  elevation  are  time  independent  due  to  a  reformulation  of  the 
WCE  by  Kinnmark  (1984)  into  a  GWCE.  Matrices  resulting  from  an  application  of  the  FE  method 
to  the  momentum  equations  remain  time  dependent,  but  this  time  dependence  is  overcome  by 
lumping  these  matrices  (that  is,  summing  off-diagonal  terms  and  adding  the  sum  to  the  diagonal 
term),  which  causes  only  a  slight  degradation  in  accuracy.  As  a  consequence,  velocities  are  computed 
by  solving  a  trivial  tridiagonal  system  of  equations  (Lynch  and  Gray  1979). 

One  final  consideration  in  selecting  a  numerical  solution  algorithm  is  grid  flexibility.  To  date, 
the  finite  difference  method  is  not  readily  amenable  to  providing  high  levels  of  refinement  in 
localized  areas.  Restrictions  on  acceptable  grid  skewness  and  maximum  cell-to-cell  size  ratios  (i.e., 
Heath  et  al.  1990;  Celia  and  Gray  1992)  often  result  in  finite  difference  grids  that  are  over-refined, 
causing  unnecessary  computational  burden.  To  represent  detailed  coastline  geometry,  coordinate 
transformation  techniques  were  developed  in  the  context  of  finite  difference  schemes,  but  these 
methods  still  do  not  allow  for  increased  resolution  in  localized  areas. 

In  contrast,  the  advantage  of  the  FE  approach  lies  in  its  tremendous  flexibility  of  nodal  placement 
to  represent  the  complexity  of  the  shoreline  and  provide  varying  degrees  of  resolution  throughout 
the  model  domain  as  warranted.  In  particular,  FE  schemes  based  on  triangular  elements  result  in 
optimal  flexibility  for  achieving  local  refinement  (e.g.,  Luettich  et  al.  1992;  Westerink  et  al.  1995; 
Blain  et  al.  1998).  With  the  ever  increasing  size  and  complexity  of  shallow-water  problems,  the 
degree  of  mesh  flexibility  inherent  in  the  discretization  strategy  directly  affects  the  efficiency  and 
accuracy  of  the  hydrodynamic  model. 


2.3  The  ADCIRC-2DDI  Hydrodynamic  Model 

A  careful  evaluation  of  the  accuracy,  efficiency,  and  grid  flexibility  characteristics  of 
current  numerical  solution  algorithms  for  the  shallow-water  equations  has  led  to  selection  of  a 
GWCE  FE-based  formulation  of  the  shallow-water  equations.  One  such  hydrodynamic  model  is 
ADCIRC-2DDI,  the  depth-integrated  portion  of  a  system  of  2D  and  three-dimensional  (3D)  codes 
named  ADCIRC  (ADvanced  CIRCulation  Model  for  Shelves,  Coasts,  and  Estuaries),  developed  by 
Luettich  et  al.  (1992)  and  Westerink  et  al.  (1992b). 

The  efficiency  and  accuracy  of  ADCIRC-2DDI  are  well  understood  due  to  extensive  numerical 
testing  and  analysis  of  the  model  code.  The  algorithms  that  comprise  ADCIRC-2DDI  effectively 
minimize  the  required  number  of  degrees  of  freedom  for  a  desired  level  of  accuracy,  show  good 
stability  characteristics,  generate  no  spurious  artificial  modes,  utilize  minimum  inherent  artificial 
numerical  damping,  efficiently  separate  the  partial  differential  equations  into  systems  of  algebraic 
equations  with  time-independent  matrices,  and  are  capable  of  running  months  to  years  of  simulation 
while  providing  detailed  computations  of  the  circulation  patterns  within  a  water  body  (Luettich  et 
al.  1992).  Other  studies  that  address  specific  characteristics  of  the  ADCIRC-2DDI  model  formulation 
are  the  investigation  of  mass  conservation  properties  of  the  GWCE  formulation  by  Kolar  et  al. 
(1994b)  and  the  accuracy  of  the  boundary  condition  formulation  within  ADCIRC-2DDI  (Westerink 
et  al.  1994b). 

The  ADCIRC-2DDI  hydrodynamic  model  produces  computations  that  are  in  agreement  with 
two  established  benchmark  problems,  the  quarter  annular  test  case  (Gray  1982)  and  the  North  Sea/ 
English  Channel  system  (Werner  and  Lynch  1988).  The  ADCIRC-2DDI  model  has  also  been  applied 
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to  a  number  of  field  studies  with  excellent  results  (Westerink  et  al.  1989,  1992  a,b,  1994a;  Kolar 
1994  a,b).  Detailed  documentation  (Luettich  et  al.  1992),  as  well  as  an  instructive  commentary 
regarding  appropriate  setup  of  the  ADCIRC-2DDI  model  (Westerink  et  al.  1992b),  facilitates  its 
application  in  the  coastal  ocean. 


2.3.1  Governing  Equations 

The  generalized  wave-continuity  equation  formulation  in  ADCIRC-2DDI  is  based  on  the 
well-known,  shallow-water  equations  (Le  Mehaute  1976;  Kinnmark  1984).  The  primitive  form  of 
the  shallow-water  equations  is  derived  by  averaging  conservation  laws  of  mass  and  momentum 
over  the  timescale  of  turbulent  fluctuations  and  ocean  depth.  Turbulent  fluctuations  at  the  microscale 
present  in  all  hydrodynamic  flows  are  characterized  by  spatial  and  temporal  variations  in  velocity 
and  pressure  fields  and  are  a  mechanism  for  momentum  transfer.  Time  averaging  the  conservation 
laws  generates  macroscale  quantities  representing  microscale  fluctuations  that  must  be  parameterized. 
The  scale  for  time  averaging  is  selected  to  be  long  enough  to  span  a  statistically  significant  sample 
of  turbulent  fluctuations,  and  yet  short  enough  so  that  macroscopic  variations  of  the  averaged 
quantity  are  not  included  (Gray  et  al.  1993).  Within  the  shallow-water  equations,  turbulent  fluctuations 
are  characterized  by  momentum  diffusion  terms. 


Depth  averaging  the  conservation  laws  reduces  a  3D  problem  to  2D,  leaving  as  unknowns  the 
surface  water  elevation  t,  and  the  depth-averaged  lateral  velocities  U  and  V.  For  flows  having  small, 
vertical  velocity  gradients  (e.g.,  well-mixed  systems)  or  environments  where  lateral  flows  are  quite 
large  in  comparison  to  vertical  velocities  (e.g.,  nearly  horizontal  flow),  application  of  the  shallow- 
water  equations  is  accepted  practice. 

Derivation  of  the  shallow-water  equations  from  the  time-  and  depth-averaged  conservation  laws 
involves  the  Boussinesq  approximation  in  addition  to  hydrostatic  and  incompressibility  assumptions. 

In  practice,  incompressibility  of  a  fluid  implies  that  density  variations  seen  while  moving  with  the 

Dp 

fluid  are  negligible  (i.e.,  =  0  ).  Furthermore,  the  Boussinesq  approximation  specifies  density  as 

a  constant  value,  p  =  p0,  except  when  gradients  of  density  are  considered.  In  flows  appropriate  for 
2D,  depth-averaged  equations,  vertical  acceleration  is  assumed  negligible.  Consequently,  the  momentum 
equation  over  the  vertical  reduces  to  a  balance  between  pressure  and  gravitational  forces.  This  is 
commonly  referred  to  as  the  hydrostatic  assumption.  Other  simplifications  implicit  within  the 
shallow-water  equations  include  neglecting  changes  in  the  position  of  the  ocean  floor  with  respect 
to  time  and  excluding  mass  exchanges  with  the  environment  (i.e.,  evaporation,  precipitation,  overland 
flow,  and  groundwater  interactions). 


For  a  Cartesian  coordinate  system,  the  conservative  form  of  the  shallow-water  equations  are  written: 


dUH  dVH 
dt  +  dx  +  dy 


(2.1) 


dUH  dUUH  dUVH  ,TFrr  d 

+  - +  — - fVH  =  -H— 

dt  dx  dy  dx 


Ps 


+  g(t-ari) 
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dVH  dVUH  dWH  rTTTT  TTd 
— T  +  — —  +  — — +fUH  =  -H— 
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where  t  represents  time,  x,  y  are  the  Cartesian  coordinate  directions,  t.  is  the  free  surface  elevation 
relative  to  the  geoid,  U,  V  are  the  depth-averaged  horizontal  velocities,  H  =  £  +  h  is  the  total  water 
column  depth,  h  is  the  bathymetric  depth  relative  to  the  geoid,  /  is  the  Coriolis  parameter,  ps  is  the 
atmospheric  pressure  at  the  free  surface,  g  is  the  acceleration  due  to  gravity,  a  is  the  Earth  elasticity 
factor,  ti  is  the  Newtonian  equilibrium  tide  potential,  p0  is  the  reference  density  of  water,  Mx,  My 
represents  the  depth-integrated  horizontal  momentum  diffusion,  Dx,  Dy  are  the  depth-integrated 
horizontal  momentum  dispersion  terms,  Bx,  By  are  the  depth-integrated  baroclinic  forcings,  and  xsx, 
xSy  are  the  applied  free  surface  stresses.  Further  justification  regarding  the  appropriateness  of  these 
equations  in  modeling  tidal  and  atmospheric  forces  flows  is  provided  by  Blumberg  and  Mellor 
(1987),  Westerink  et  al.  (1989),  and  Luettich  et  al.  (1992). 

In  operator  notation,  a  relationship  between  the  conservative  and  nonconservative  momentum 
equations  is  given  (Kolar  1992): 


MNC  =  ±(MC-vC),  (2.4) 

where  are  the  nonconservative  momentum  equations,  represents  the  conservative  momentum 
Eqs.  (2.2)  and  (2.3),  v  is  the  horizontal  velocity  vector,  and  C  is  the  primitive  continuity  Eq.  (2.1). 
Substituting  Eqs.  (2.1)-(2.3)  into  Eq.  (2.4)  leads  to  a  reformulation  of  the  momentum  equations  into 
primitive,  nonconservative  form: 


dt  dx  dy 


d 

Ps  , , 

1 

dx 

— +  g(^-ari) 
Vo 

+  H 
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dt  dx  dy 
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(2.6) 


A  rigorous  derivation  of  Eqs.  (2.1)-(2.6)  is  presented  by  Kolar  (1992)  and  will  not  be  repeated  here. 

Implementation  of  a  standard  quadratic  parameterization  for  bottom  stress  and  the  neglect 
of  baroclinic  and  lateral  diffusion/dispersion  terms  leads  to  a  modified  form  of  the  primitive, 
nonconservative,  shallow-water  equations: 


dUH  dVH 
dt  dx  dy 


(2.7) 
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with  x*  given  by  the  expression  Cf  (U2  +  V2)y2/H,  for  Coequal  to  the  bottom  friction  coefficient. 
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Recall  that  one  objective  of  the  ADCIRC-2DDI  hydrodynamic  model  is  to  be  able  to  consider 
very  large  domain  problems.  The  use  of  large  domains  coupled  with  a  recognition  that  tidal  forcing 
is  a  global  phenomena  necessitates  the  inclusion  of  effects  caused  by  curvature  of  the  Earth’s 
surface.  Thus,  shallow-water  Eqs.  (2.7)-(2.9)  are  recast  into  spherical  coordinates  (Flather  1988; 
Kolar  et  al.  1994a): 


dt,  1  dukH  cos  4> 

dt  +  R  cos  <)>  dX  +  d<|)  "  ^ 
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(2-12) 


where  X,  <J>  are  degrees  longitude  (east  of  Greenwich  positive)  and  degrees  latitude  (north  of  the 
equator  positive),  u^,  are  depth-averaged  velocities  in  spherical  coordinates,  R  is  the  radius  of 
the  Earth,  /  is  given  by  2Q  sin  <j>,  Q  is  the  angular  speed  of  the  Earth,  and  xs^,  are  the  applied 
free  surface  stresses. 

A  practical  expression  for  the  Newtonian  equilibrium  tide  potential  is  given  by  Reid  (1990)  as: 


X](k,  <M)  =  2  cjn  fjn  (f0)  Ty(4»)  cos  [2n(t-t0)/Tjn  +jX  +  vjn(t0)] , 


(2.13) 


where  t  is  time  relative  to  the  reference  time,  to,  Cjn  is  a  constant  characterizing  the  amplitude  of 
tidal  constituent  n  of  species  j,  fjn  is  a  time-dependent  nodal  factor,  Vjn  is  the  time-dependent 
astronomical  argument,  j  =  0,  1,  2  are  the  tidal  species  (  j  =  0  declinational;  j  =  1  diurnal,  j  =  2 
semidiurnal),  Lq  =  3sin2<|>  - 1,  L\  =  sin(2<J>),  Z,2  =  cos2(<|>),  and  Tjn  is  the  period  of  constituent  n  for 
species  j.  Values  for  C-jn  are  presented  by  Reid  (1990)  and  the  Earth  elasticity  factor  is  often  taken 
as  0.69  for  all  tidal  constituents  (Schwiderski  1980;  Hendershott  1981),  although  its  value  has  been 
shown  to  be  slightly  constituent-dependent,  ranging  between  0.693  and  0.736  (Wahr  1981;  Woodworth 
1990). 

The  Earth’s  curvature  must  be  accounted  for,  not  only  in  the  governing  equations,  but  also  in 
the  FE  discretization  (Kolar  et  al.  1994a).  In  the  FE  method,  the  solution  to  the  governing  equations 
is  approximated  using  interpolating  functions.  These  interpolating  functions  are  defined  over 
elements  and  these  elements  are  most  often  cast  within  the  framework  of  a  Cartesian  coordinate 
system.  To  conveniently  implement  the  FE  method,  governing  equations  in  spherical  coordinates 
are  projected  onto  a  planar  surface  using  cartographic  projection  techniques.  Mapping  spherical 
Eqs.  (2.10)-(2.12)  to  a  rectilinear  coordinate  system  is  accomplished  using  a  Carte  Parallelogramatique 
Projection  (CP)  (Pearson  1990): 


x'  =  R(k  -  Xa)  cos  <j>0 


(2.14) 


(2.15) 


Slain  and  Rogers 


where  (X0,  <t>o)  is  the  center  point  of  the  projection.  An  application  of  the  CP  projection  to 
Eqs.  (2.10)-(2.12)  yields  shallow-water  equations  in  primitive,  nonconservative  form  expressed  in 
the  CP  coordinate  system: 
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As  discussed  in  Sec.  2.2,  utilizing  the  FE  method  to  resolve  spatial  dependencies  in  the  primitive 
shallow-water  equations  leads  to  inaccurate  solutions  with  severe  artificial,  near  2  •  Ax  modes  (Gray 
1982;  Lynch  1983;  Westerink  et  al.  1989).  Reformulation  of  the  primitive  equations  into  a  GWCE 
gives  highly  accurate,  noise-free,  FE-based  solutions  to  the  shallow-water  equations  (Lynch  and 
Gray  1979;  Kinnmark  1984).  The  high  accuracy  of  GWCE-based  FE  solutions  is  a  result  of  their 
excellent  numerical  amplitude  and  phase  propagation  characteristics.  In  fact,  Fourier  analysis 
indicates  that  for  constant-depth  water  using  linear  interpolation,  a  linear  tidal  wave  resolved  with 
25  nodes  per  wavelength  is  more  than  adequately  represented  over  the  range  of  Courant  numbers, 
C  =  \fgh  At/Ax,  less  than  or  equal  to  one  (Luettich  et  al.  1992).  Furthermore,  the  monotonic  dispersion 
behavior  of  GWCE-based  FE  solutions  avoids  generating  artificial,  near  2  •  Ax  modes  that  have 
plagued  the  primitive  FE-based  solutions  on  the  interior  and  at  the  boundary  (Platzman  1981; 
Foreman  1983;  Westerink  et  al.  1994b).  Note  that  the  monotonic  dispersion  behavior  of  GWCE- 
based  FE  solutions  is  very  similar  to  that  associated  with  staggered  finite  difference  solutions  to 
the  primitive  shallow-water  equations  (Westerink  and  Gray  1991).  GWCE-based  FE  solutions  to  the 
shallow-water  equations  allow  extremely  flexible  spatial  discretizations  that  effectively  minimize 
the  discrete  size  of  any  problem  (LeProvost  and  Vincent  1986;  Foreman  1988;  Vincent  and  LeProvost 
1988;  Westerink  et  al.  1994a,  1995;  Luettich  and  Westerink  1995). 


Derivation  of  the  GWCE  is  presented  concisely  using  the  operator  notation  invoked  by  Kinnmark 
(1984)  and  Kolar  (1992): 


GWCE  = 


dC 

dt 


+  t0C- 


V  •  Mc, 


(2.19) 


where  GWCE  is  the  generalized  wave-continuity  equation  and  xa  is  a  nonphysical  constant  in  time 
and  space  that  controls  the  balance  between  primitive  and  wave  equation  formulations  (Lynch  and 
Gray  1979;  Kinnmark  1984;  Luettich  et  al.  1992;  Kolar  et  al.  1994b).  Substituting  continuity 
Eq.  (2.16)  and  the  conservative  form  of  momentum  Eqs.  (2.17)  and  (2.18)  into  Eq.  (2.19),  the 
GWCE  in  the  CP  coordinate  system  is: 
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Advective  terms  in  the  GWCE,  written  in  nonconservative  form,  improve  global  and  local  mass 
conservation  (Kolar  et  al.  1994b),  as  well  as  numerical  stability,  especially  for  advection  dominant 
flows  (Westerink  et  al.  1992b).  The  ADCIRC-2DDI  hydrodynamic  model  solves  the  GWCE,  Eq.  (2.20), 
in  conjunction  with  the  primitive  momentum  equations  in  nonconservative  form,  Eqs.  (2.17)  and  (2.18). 


2.5.2  Numerical  Solution 

Numerical  discretization  of  GWCE  Eq.  (2.20)  and  momentum  Eqs.  (2.17)  and  (2.18)  is  implemented 
in  three  stages.  First,  the  symmetrical  weak  weighted  residual  (SWWR)  statements  for  Eqs.  (2.17), 
(2.18),  and  (2.20)  are  developed.  This  procedure  is  based  on  a  standard  Galerkin  FE  formulation 
(Becker  et  al.  1981;  Celia  and  Gray  1992).  One  consequence  of  the  SWWR  is  that  the  order  of  the 
derivatives  in  the  governing  equations  is  reduced,  leading  to  a  requirement  of  only  C0  functional 
continuity  (i.e.,  only  the  interpolating  functions  themselves  and  not  their  derivatives  need  be  continuous 
between  discrete  points).  Next,  the  stable  and  accurate  time  discretization  strategies  of  Kinnmark 
(1984)  and  Werner  and  Lynch  (1989)  are  implemented.  A  variably  weighted  three-time-level  implicit 
scheme  is  applied  to  linear  terms  in  the  GWCE.  The  nonlinear  Coriolis,  atmospheric  pressure,  and 
tidal  potential  terms  are  all  treated  explicitly.  Alternatively,  advective  terms  within  the  GWCE  are 
evaluated  at  two  known  time  levels.  This  time  discretization  results  in  a  system  of  linear 
algebraic  equations  associated  with  the  GWCE,  which  is  solved  for  unknown  elevations.  In  the 
momentum  equations,  a  Crank-Nicolson  two-time-level  implicit  scheme  is  applied  to  all  terms 
except  the  bottom  friction  and  advective  terms,  which  are  treated  explicitly. 

The  final  step  in  the  numerical  discretization  scheme  is  approximation  of  the  spatial  domain 
using  the  FE  method.  Variables  are  expanded  using  a  C0  interpolation  basis  over  three-node,  linear 
triangular  elements.  Elemental  equations  are  summed  over  the  global  domain  and  inter-element  Ca 
functional  continuity  is  enforced.  Details  of  the  FE  implementation  on  a  term-by-term  basis  is 
presented  by  Luettich  et  al.  (1992)  and  Westerink  et  al.  (1992b). 

The  fully  discretized  model  equations  are  written  in  matrix  notation  (Luettich  et  al.  1992): 


mGWCE  %k  +  1  _  pGWCE 


(2.21) 
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(2.22) 


(2.23) 


where  mgwce  is  the  banded,  time-independent  mass  matrix  in  the  GWCE  equation,  %k  +  1  is  the 
vector  of  unknown  surface  elevations  at  time  level  k  +  1,  pGWCE  is  the  load  vector  of  known 
forcings  in  the  GWCE  equation,  Mxxme,  Mx$me,  M$xme,  and  M^me  are  the  time-dependent, 
lumped  mass  matrices  in  the  X,<f>  directions,  Uk+1,  Vk  +  1  are  the  vectors  of  unknown  velocity 
components  at  the  k  +  1  time  level,  and  p^ME,  p§ME  are  ^jje  ]0ad  vectors  of  known  forcings  for  the 
momentum  equations. 


Elevation  boundary  conditions  are  enforced  within  the  load  vector  pGWCE  of  the  GWCE  Eq.  (2.21) 
and  zero  normal  velocity  boundary  conditions  are  enforced  in  the  momentum  Eqs.  (2.22)  and 
(2.23).  Westerink  et  al.  (1994b)  have  shown  that  solutions  to  the  GWCE  equation  are  insensitive 
to  this  standard  boundary  condition  formulation. 


The  decoupled  discrete  GWCE  and  momentum  Eqs.  (2.21)-(2.23)  lead  to  a  sequential  solution 
procedure.  GWCE  Eq.  (2.21)  is  solved  at  each  timestep  for  the  surface  water  elevations  t,k  +  *.  The 
GWCE  mass  matrix  MGWCE  is  time-independent,  so  it  is  assembled  and  decomposed  only  once. 
The  banded  structure  of  this  matrix  is  not  utilized  by  the  iterative  preconditioned  conjugate  gradient 
matrix  solver  used  for  computations  performed  in  this  work  (Press  et  al.  1986;  Kincaid  and  Cheney 
1991).  A  preconditioned  conjugate  gradient  solver  is  implemented  because  of  its  efficiency  at 
minimizing  memory  requirements  within  ADCIRC-2DDI  for  larger  problems.  The  load  vector  PGWCE 
is  updated  at  each  timestep  with  newly  computed  surface  water  elevations  and  velocities  from  the 
previous  timestep. 

After  solving  the  GWCE  for  surface  water  elevations,  these  computed  elevations  are  substituted 
into  momentum  Eqs.  (2.22)  and  (2.23)  prior  to  solution  for  the  velocity  components  and  Uk  +  1  and 
V*  +  1.  The  time-dependent  mass  matrices  in  the  momentum  equations  Mxxme,  mx^me,  m^xme, 
and  M^me  are  lumped  to  yield  diagonal  matrices  that  require  trivial  solution.  The  lumping  procedure 
applied  here  does  not  introduce  significant  errors  as  shown  by  Lynch  and  Gray  (1979). 

The  numerical  solution  algorithm  just  described  for  the  ADCIRC-2DDI  model  is  implemented 
in  fully  vectorized  form.  A  consequence  of  this  solution  procedure  is  a  highly  efficient  code  in 
terms  of  central  processing  unit  (CPU)  requirements  per  node.  This  efficiency  is  largely  due  to  the 
fact  that  GWCE-based  FE  solutions  to  the  shallow-water  equations  allow  for  extremely  flexible 
spatial  discretizations  that  result  in  minimization  of  the  discrete  size  of  any  problem  (LeProvost  and 
Vincent  1986;  Foreman  1988;  Vincent  and  LeProvost  1988;  Westerink  et  al.  1992  a,b,  1994a,  1995). 


2.3.3  Summary 

Numerical  formulation  of  the  hydrodynamic  model  is  critical  to  achieving  successful  predictions 
of  tidal  circulation  on  the  continental  shelf.  The  depth-averaged  shallow-water  equations  are  an 
appropriate  theoretical  framework  for  the  barotropic  hydrodynamic  model  with  additional  requirements 
of  accuracy,  efficiency,  and  grid  flexibility  imposed  on  the  numerical  solution  technique.  In  considering 
available  numerical  solution  techniques  for  the  shallow-water  equations,  a  GWCE-based  FE  approach 
is  thought  to  be  an  optimal  solution  strategy. 
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The  hydrodynamic  model,  ADCIRC-2DDI,  solves  the  GWCE-based  FE  equations  and  has 
demonstrated  the  desired  characteristics  of  accuracy,  efficiency,  and  a  high  degree  of  grid  flexibility 
in  numerous  idealized  tests,  analyses,  and  field  applications.  For  all  tidal  simulations  presented  in 
this  work,  the  ADCIRC-2DDI  hydrodynamic  model  is  utilized.  Modularity  of  the  model  allows  the 
user  to  render  various  terms  and/or  forcing  mechanisms  within  the  model  active  or  inactive.  This 
modularity  is  found  to  be  a  necessary  asset  both  for  the  study  described  herein  and  for  future 
diagnostic  examinations  of  coastal  tidal  dynamics. 


3.0  THE  TIDAL  FLOW  FORUM  BENCHMARK  CASE 
3.1  Description 

Tidal  Flow  Forums  1  (TFF)  (Lisbon,  Portugal  1986)  and  2  (Cambridge,  MA  1988)  provide 
benchmark  test  cases  for  the  quantitative  skill  assessment  of  coastal  tidal  hydrodynamic  models. 
The  southern  North  Sea/English  Channel  is  chosen  by  TFF  organizers  as  an  appropriate  study  area 
for  several  reasons:  (1)  the  hydrodynamics  of  the  region  have  been  studied  extensively,  (2)  a  fairly 
complete  set  of  model  input  data  is  readily  available,  and  (3)  field  observations  of  both  elevation 
and  currents  are  available.  Characteristic  nonlinearity  and  spatial  variability  of  tides  in  the  region 
provide  a  good  test  for  any  tidal  hydrodynamic  model.  Furthermore,  the  small  size  of  the  domain 
makes  a  tractable  problem  for  a  variety  of  numerical  model  formulations  and  computer  resources, 
facilitating  model-model  comparisons.  Essential  information  for  model  setup  and  validation  data 
provided  by  Verboom  and  LeProvost  (1986)  include  bathymetry,  computational  meshes,  domain 
boundary  tidal  forcing,  sea  surface  height,  and  current  meter  field  data.  The  TFF  standardizes  a 
single  tide  model  application  for  the  purpose  of  validation  and  inter-comparison  of  numerical  tidal 
models. 

In  the  spirit  of  the  TFF,  model-to-model  and  model-to-data  comparisons  using  both  finite 
difference  (Ozer  and  Jamart  1988;  Jamart  and  Ozer  1989;  Yu  et  al.  1989;  Praagman  et  al.  1989  a,b) 
and  FE  (Werner  and  Lynch  1987;  Gray  et  al.  1987;  Walters  1987;  Lynch  and  Werner  1988;  Werner  and 
Lynch  1989;  Gray  1989;  Baptista  et  al.  1989;  Walters  and  Werner  1989;  Werner  1995)  models  were 
conducted.  The  TFF  study  area  has  been  modeled  previously  using  ADCIRC-2DDI  by  Luettich 
et  al.  (1992).  For  this  application  of  ADCIRC-2DDI  in  the  southern  North  Sea/English  Channel, 
Luettich  et  al.  (1992)  publish  only  a  cursory  comparison  to  TFF  data  from  a  limited  number  of 
model  simulations.  The  present  study  is  a  considerable  extension  of  that  initial  validation  effort. 
The  model  experiments  presented  herein  are  intended  to  provide  greater  insight  into  simulating  tidal 
physics  with  the  ADCIRC-2DDI  hydrodynamic  model. 


3.1.1  Model  Domain 

The  model  domain,  shown  in  Fig.  1,  extends  from  the  Dutch  and  English  coasts  in  the  southern 
North  Sea  to  the  French  and  English  coasts  at  the  western  end  of  the  English  Channel.  The  original 
bathymetric  data  provided  by  Verboom  and  LeProvost  (1986)  has  since  been  revised  to  resolve 
inconsistencies  in  the  data  (Werner  and  Lynch  1988).  The  FE  grid,  shown  in  Fig.  2,  consists  of 
911  nodes  and  1613  elements.  Nodal  spacings  for  this  mesh  (labeled  G9  indicating  a  maximum 
resolution  of  9  km)  range  between  9-15  km.  Islands  are  not  represented  in  this  bathymetric  grid; 
one  notable  omission  is  the  Isle  of  Wight  near  Christchurch. 
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Fig.  1  —  Southern  North  Sea/English  Channel  bathymetry 


3.1.2  Open  Boundary  Forcing 

Tidal  elevation  data  in  the  form  of  tidal  phases  and  amplitudes  are  provided  at  23  locations  on 
two  open-ocean  boundaries  (identified  in  Fig.  1)  for  11  constituents.  The  11  tidal  constituents 
include  six  primary  components — two  diurnal  (Oj,  K{),  four  semidiurnal  (M2,  S2,  N2,  K2),  and  five 
nonlinear  constituents  (A/4.  MS 4,  MN 4,  Af6,  2A/5g).  These  11  components  will  be  referred  to  as 
“TFF  constituents.”  This  set  of  data  is  derived  from  the  harmonic  decomposition  of  measured 
elevations  and  are  supplemented  by  the  physical  model  study  of  d’Hieres  and  LeProvost  (1979). 
The  TFF  constituent  data  are  linearly  interpolated  onto  boundary  nodes  of  the  FE  mesh  and  constitute 
the  sole  means  of  forcing  in  this  study  unless  otherwise  indicated. 

3.1.3  Bottom  Friction  Parameterization 

For  purposes  of  standardization,  TFF  organizers  suggest  a  value  for  the  Chezy  bottom  friction 
coefficient  equal  to  65  m0-5/s  (Verboom  and  LeProvost  1986).  For  a  quadratic  bottom 
friction  representation,  the  Chezy  coefficient  translates  to  a  value  of  Cy  equal  to  2.322  x  10-3.  This 
is  the  frictional  coefficient  used  in  ADCIRC-2DDI  for  those  simulations  that  utilize  a  nonlinear 
representation  of  bottom  friction. 
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3.2  Field  Observations 

Amplitudes  and  phases  for  11  tidal  constituents  constitute  the  tidal  elevation  data  provided  at 
11  coastal  locations  shown  in  Fig.  1.  These  same  TFF  tidal  constituents  are  used  for  boundary 
forcing  and  are  thought  to  be  of  high  quality  since  they  are  derived  from  long-term  observations. 
However,  as  will  be  discussed  in  Sec.  5.1.2,  these  11  tidal  constituents  do  not  provide  a  complete 
representation  of  the  true  tidal  signal;  other  constituents  are  shown  to  make  significant  contributions 
in  most  areas  of  the  domain.  The  field  data  available  to  the  TFF  organizers  was  not  adequate  to 
allow  inclusion  of  these  other  important  tidal  constituents  in  the  harmonic  analyses.  As  such,  the 
time  series  of  elevation  at  each  elevation  gauge  location  is  constructed  from  available  tidal  constituent 
information  using  the  well-known  harmonic  relation  for  tides  (Schureman  1958): 

n 

h=  2  fi  cos  (o)jt  -  Gt  +  Ei).  (3.1) 

i  =  l 

Here,  h  is  the  sea  surface  elevation,  n  is  the  number  of  tidal  constituents  included  in  the  time  series, 
/  is  the  constituent  nodal  factor,  A  is  the  constituent  amplitude  at  that  location,  co  is  the  constituent 
frequency,  t  is  time,  G  is  the  constituent  phase  at  that  location,  and  E  is  the  constituent  equilibrium 
argument  phase  correction.  Henceforth,  time  series  reconstructed  from  the  data  are  referred  to  as 
the  observed  elevation  time  series. 
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At  eight  offshore  locations  also  shown  in  Fig.  1,  velocity  data  taken  from  JONSDAP’76  (Ramster 
1977)  are  based  on  direct  measurement  at  a  single  location  in  the  water  column.  From  comparisons 
of  the  velocity  data  to  ADCIRC  model  solutions,  as  well  as  from  previous  TFF  discussions,  overall 
speed  observations  are  of  questionable  value.  Further  discussion  pertaining  to  the  velocity  data  is 
found  in  Sec.  5.2. 


4.0  MODEL  VALIDATION  EXPERIMENTS 
4.1  Description 

In  general,  model  experiments  performed  for  this  study  are  classified  according  to  two  broad 
categories  based  on  length  of  the  simulation.  In  keeping  with  the  TFF  benchmark  case,  a  preliminary 
simulation  of  3  d  from  15-17  Mar  1976  is  conducted.  The  first  2  d  serve  as  a  ramp-up  period  with 
analysis  during  the  last  24  h.  In  addition,  45-d  experiments  using  a  30-d  spin-up  period  allow  time 
for  full  dissipation  of  the  free  Helmholtz  modes  (Westerink  et  al.  1994a).  The  45-d  experiments 
include  a  15-d  period  at  the  beginning  of  the  simulation  during  which  ADCIRC’s  hyperbolic  tangent 
ramp  function  is  applied  to  the  model  forcing  (Luettich  et  al.  1992).  Application  of  the  hyperbolic 
tangent  function  results  in  zero  forcing  at  the  initiation  of  the  experiment  ( t  =  0)  and  full-strength 
forcing  at  the  end  of  the  15-d  ramp-up  period  (r  =  15  d).  Some  minor  differences  are  observed 
between  the  45-d  and  3-d  experiments,  so  a  45-d  duration  is  used  for  all  short-term  model  experiments 
presented.  The  model  reference  time,  used  to  compute  the  tidal  nodal  factors  and  equilibrium 
arguments,  is  defined  42  d  into  the  short  duration  simulations;  the  date  corresponding  to  this 
reference  time  is  15  Mar  1976.  Model  solutions  from  these  short  simulations  are  compared  to 
observed  time  series  of  elevation  and  velocity. 

Harmonic  decomposition  based  on  the  least  squares  analysis  method  (e.g.,  van  Ette  and  Schoemaker 
1967;  Foreman  1977)  is  implemented  to  decompose  the  model-computed  time  series  into  tidal 
constituent  frequencies.  Harmonic  analyses  of  model  solutions  for  short  (45-d)  simulations  are 
unsuccessful,  largely  because  of  the  inability  of  the  method  to  distinguish  the  K-2  and  52  frequencies 
using  an  abbreviated  time  signal.  In  fact,  Godin  (1972)  recommends  a  time  series  length  greater 
than  183  d  to  accurately  extract  the  K2  and  52  tidal  constituents.  Therefore,  a  set  of  long-term 
simulations  of  205  d  in  length  was  designed.  These  model  experiments  have  15  d  for  an  initial 
ramp-up  period  and  the  remaining  190  d  of  the  simulation  are  used  in  the  harmonic  analyses  for 
56  primary,  compound,  and  overtide  constituents.  The  harmonically  decomposed  time  series  of 
elevations  are  then  utilized  in  model-to-model  and  model-to-data  comparisons  for  individual  tidal 
constituents  or  groups  of  tidal  constituents. 


4.2  Quantitative  Analyses 

4.2.1  Error  Measures 

In  analyzing  model  performance,  two  error  measures  are  used  for  the  time  series  comparisons:  the 
root-mean-square  (RMS)  error,  erms,  and  a  proportional  error  e2.  The  RMS  error  is  expressed: 


erms 


^  ixmod  xobs ) 

j=  1 


(4.1) 
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and  the  proportional  error  is  computed 


^  ( xmod  xobs) 


e2  =  ^i- 


obs 


;  =  i 


(4.2) 


Here,  xmod  and  x0ys  are  the  modeled  and  observed  elevation,  current  speed,  or  direction,  respectively, 
and  n  is  the  number  of  observation  stations.  Computations  of  these  two  errors  are  based  on  comparisons 
between  model-computed  elevation,  speed,  direction,  and  observations  of  the  same  fields  over  the 
last  3  d  of  the  model  simulation,  15-17  Mar  1976.  These  error  measures  are  intended  to  gauge 
accuracy  of  the  ADCIRC-2DDI  tidal  predictions. 

Errors  for  individual  tidal  constituents  are  quantified  by  distance  in  the  complex  plane  D, 
calculated  (Foreman  and  Henry  1993): 


D  —  \J (A 0 cos PQ  A  m cosPJ  +  (A o Pq  Aflj sin Pm )  ,  (4.3) 


where  Aa,  Am,  Pa,  and  Pm  are  the  observed  and  modeled  amplitudes  and  phases,  respectively. 


4.2.2  Mass  Balance  Check 

An  important  property  of  any  coastal  circulation  model  is  mass  conservation.  For  several 
model  experiments,  global  mass  balance  checks  are  conducted  by  comparing  changes  in  water 
volume  within  the  domain  to  the  flux  in/out  of  the  domain.  The  methodology  used  for  these 
calculations  is  described  in  detail  by  Kolar  et  al.  (1994b). 


4.3  Base  Simulation 

Taking  into  consideration  the  suggestions  of  TFF  organizers  and  participants,  as  well  as  using 
experience  gained  from  numerous  ADCIRC  model  simulations,  a  base  simulation  is  created  that 
offers  an  ideal  combination  of  simplicity  and  performance  (accuracy)  without  tuning  model  parameters. 
The  model  equations  implemented  are  fully  nonlinear  and  include  a  quadratic  parameterization  of 
nonlinear  bottom  friction  forces,  nonlinear  convective  acceleration  terms  (both  space  and  time 
derivatives),  and  nonlinear  finite  amplitude  terms.  The  quadratic  nonlinear  friction  factor  Cy  is  set 
equal  to  2.322  x  10-3  and  the  GWCE  weighting  factor  is  defined  by  xa  =  0.0006.  The  meaning  for 
these  parameters  was  given  in  Sec.  2.3.1.  Along  the  open  boundary,  11  tidal  constituents  (six 
primary,  five  nonlinear)  are  forced.  For  all  simulations,  no  normal  flow,  free  tangential  slip  land 
boundaries  are  used.  Tidal  potential  forcing  is  not  activated  and  the  Coriolis  parameter  is  spatially 
constant,  a  consequence  of  the  Cartesian  coordinate  system  that  is  employed.  Wind  forcing  is  not 
included  in  this  study. 

The  coefficient  for  lateral  eddy  viscosity  is  equal  to  zero,  and  a  minimum  depth  of  10  m  is 
specified  to  prevent  drying  of  mesh  elements.  A  complete  description  of  the  ADCIRC-2DDI  model 
equations  was  provided  previously  in  Sec.  2.3. 


16 


Blain  and  Rogers 


5.0  RESULTS  FOR  THE  BASE  SIMULATION 

Semidiurnal  frequencies  dominate  the  tidal  response  in  the  southern  North  Sea/English  Channel. 
For  the  primary  semidiurnal  tides,  which  include  A/2,  S2,  K2,  and  N2,  one  amphidrome  is  present 
in  the  southern  North  Sea  and  one  virtual  amphidrome  is  evident  in  the  English  Channel  near 
Christchurch.  Relative  to  the  southern  North  Sea,  sea  surface  heights  associated  with  the  semidiurnal 
tides  tend  to  be  higher  on  the  continental  side  of  the  water  body  and  in  the  English  Channel.  These 
characteristics  are  shown  in  Fig.  3a  and  b  for  the  dominant  M2  tide.  Diurnal  tides  tend  to  be  larger 
in  the  North  Sea  than  in  the  English  Channel  and  exhibit  no  real  amphidromes.  The  quarter-diurnal 
tides  have  amphidromes  evenly  distributed  throughout  the  model  domain  with  two  in  the  English 
Channel  and  one  or  more  in  the  southern  North  Sea. 

The  tidal  response  at  three  transects  across  the  model  domain  are  shown  in  Fig.  4a,  which 
includes  the  underlying  bathymetry  and  Fig.  4b,  which  places  the  transects  in  the  Cartesian  coordinate 
framework.  Amplitudes  of  the  dominant  diurnal,  semidiurnal,  and  M2  overtide  frequencies  across 
each  transect  are  compared  in  Fig.  5a— c.  Minima  observed  in  these  plots  are  indicative  of  proximity 
to  an  amphidrome  or  virtual  amphidrome.  In  general,  the  higher  frequency  nonlinear  tides  exhibit 
greater  spatial  variability  and  sensitivity  to  bathymetry  as  expected. 

Across  transect  1,  a  virtual  amphidrome  of  the  Kx  tide  is  apparent  near  the  channel  constriction 
at  the  Dover  Strait.  The  dominant  M2  frequency  and  its  A/4  overtide  are  both  larger  at  the  southern 
end  of  this  transect,  which  is  located  in  waters  slightly  shallower  than  those  at  the  northern  reach 
of  the  transect.  Again,  local  lows  in  the  A/g  amplitude  indicate  positions  of  two  amphidromes. 
Amplitudes  of  the  nonlinear  tides,  A/4  and  A/g,  are  approximately  10-12%  of  the  A/2  tidal  amplitude; 
these  nonlinear  tides  are  largest  near  the  Dover  Strait,  suggesting  possible  excitation  by  the  geometry 
at  the  constriction. 

Transects  2  and  3  clearly  show  that  the  dominant  astronomical  tides  are  larger  near  the  French 
side  of  the  English  Channel.  This  is  attributed  to  a  Kelvin  wave  that  propagates  northward  through 
the  channel  (Proudman  1953).  Across  transect  2,  the  A/4  and  A/g  constituents  are  excited  on  the 
shallow,  broad  shelf  along  the  northern  coast;  the  same  is  true  for  the  Af4  constituent  on  the  mildly 
sloped  northern  coast  of  transect  3.  In  contrast,  the  A/g  constituent  over  transect  3  is  small  due  to 
the  presence  of  an  amphidrome. 

No  overall  conclusions  can  be  drawn  regarding  the  generation  of  nonlinear  tides  in  the  North 
Sea/English  Channel  model  domain.  Depths  everywhere  in  the  model  domain  are  considered  shallow 
in  relation  to  tidal  dynamics  and,  thus,  nonlinear  tides  are  largely  affected  by  local  bathymetric 
changes. 


5.1  Elevation  Time  Series  Comparisons 

5.1.1  Observations  vs.  Raw  Model  Predictions 

Initially,  raw  model  computations  for  the  base  simulation  are  compared  to  observations.  In 
general,  the  observed  time  series  for  elevation  that  contains  only  11  tidal  constituents  compares 
favorably  with  raw  model  predictions  comprised  of  an  infinite  number  of  internally  generated 
constituents.  Table  1  presents  an  RMS  error  for  each  of  the  elevation  stations.  The  full  set  of  tidal 
elevation  time  series  model-data  comparisons  using  the  base  simulation,  raw  model  predictions  is 
found  in  App.  A.  Differences  evident  in  these  time  series  warrant  further  investigation  and  a  truer 
model-data  comparison  using  “filtered”  model  computations. 
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Fig.  3  —  M2  base  simulation,  (a)  amplitudes  (meters)  and  (b)  phases  (degrees) 
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Fig.  5  —  Constituent  amplitude  variations  across  transects  (a)  1,  (b)  2, 

and  (c)  3 


20 


Blain  and  Rogers 

Table  1  —  RMS  Errors  for  Raw  Elevation  Model-Data  Comparisons 


GAUGE  LOCATION 

RMS  ERROR  (cm) 

GAUGE  LOCATION 

RMS  ERROR  (cm) 

St.  Malo 

39.6 

Hoek  van  Holland 

16.6 

Cherbourg 

22.4 

Walton 

36.9 

Dieppe 

55.5 

Dover 

28.4 

Boulogne 

40.7 

Christchurch 

24.9 

Calais 

49.6 

Lowestoft 

8.4 

Zeebrugge 

34.3 

MEAN 

32.5 

5.1.2  Observations  us.  “Filtered”  Model  Predictions 

A  fair  model-data  comparison  involves  only  that  portion  of  the  model  response  due  to  the  same 
11  constituents  comprising  the  tidal  elevation  gauge  observations.  Long-term  simulations  harmonically 
decomposed  for  56  tidal  constituents  are  used  to  reconstruct  a  time  series  based  on  only  these 
11TFF  constituents.  Time  series  reconstruction  follows  from  Eq.  (3.1).  The  reconstructed  time 
series  essentially  “filters”  all  other  tidal  frequencies  from  the  model-computed  response.  Use  of  the 
filtered  model  solution  results  in  more  valid  model-data  comparisons.  Overall  agreement  between 
the  filtered  model  predictions  and  the  data  is  improved  as  seen  from  the  RMS  errors  presented  in 
Table  2.  A  complete  set  of  filtered  tidal  elevation  model-data  comparisons  at  all  1 1  stations  is  found 
in  App.  B.  For  the  period  15-18  Mar,  an  average  proportional  error,  e2  =  0.010,  is  calculated  for 
10  of  the  11  stations.  The  remaining  station,  Christchurch,  is  excluded  from  this  average  because 
of  its  obviously  deviant  proportional  error  of  0.086.  Most  likely,  the  poor  agreement  at  Christchurch 
can  be  attributed  to  the  station’s  location  near  a  semidiurnal  “virtual”  amphidrome.  Near  this  virtual 
amphidrome,  semidiurnal  tidal  amplitudes  are  small  and  large,  spatial  gradients  of  phase  are  evident 
(refer  to  Fig.  3a  and  b). 

The  time  series  for  elevation  gauges  at  Calais,  Hoek  van  Holland,  and  Christchurch  are  presented 
in  Fig.  6a-c.  Improvements  in  the  agreement  between  model  and  data  suggest  that  compared  to 
errors  previously  computed  for  the  raw  model  solution,  constituents  excluded  from  the  data  are  not 
negligible.  An  elevation  time  series  of  the  excluded  tidal  constituents  at  each  of  stations  Calais, 
Hoek  van  Holland,  and  Christchurch  (shown  in  Fig.  7a-c)  confirms  that  these  constituents  are 


Table  2  —  RMS  Errors  for  Filtered  Elevation  Model-Data  Comparisons 


GAUGE  LOCATION 

RMS  ERROR  (cm) 

GAUGE  LOCATION 

RMS  ERROR  (cm) 

St.  Malo 

19.2 

Hoek  van  Holland 

12.7 

Cherbourg 

13.9 

Walton 

26.4 

Dieppe 

29.8 

Dover 

14.0 

Boulogne 

14.9 

Christchurch 

14.4 

Calais 

30.7 

Lowestoft 

8.8 

Zeebrugge 

16.2 

MEAN 

18.3 

ELEVATION  (m)  ELEVATION  (m) 
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Fig.  6  —  Base  simulation  tidal  elevations  at  tidal  gauges  (a)  Calais,  (b)  Hoek  van  Holland,  and 

(c)  Christchurch 
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indeed  significant.  From  the  harmonic  analysis,  constituents  M2  ancl  V-2  (not  included  in  the  TFF 
data)  are  typically  larger  than  the  smallest  signals  associated  with  the  11  TFF  constituents.  Other 
investigators  have  made  similar  observations  based  on  computations  from  their  models  (e.g.,  Werner 
and  Lynch  1989;  Jamart  and  Ozer  1989).  Even  though  errors  with  respect  to  the  measured  data  are 
smaller  for  filtered  model  results,  one  should  keep  in  mind  that  in  all  likelihood,  the  raw  model  time 
series  predictions  are  closer  to  reality.  True  sea  surface  elevations  are  unfiltered  and  contain  contributions 
from  a  wide  band  of  frequencies.  For  this  reason,  raw  model  solutions  are  used  for  the  sensitivity 
analyses  that  follow  in  Sec.  6.0. 

Apart  from  the  elevation  time  series  comparisons,  the  11  TFF  constituents  are  compared  individually 
to  the  observations.  Appendix  C  contains  figures  that  depict,  for  each  station  and  tidal  constituent, 
the  amplitude  and  phase  errors  and  a  calculated  value  for  the  error  in  the  complex  plane  D.  Con¬ 
stituent  phase  errors,  on  average,  are  larger  than  error  associated  with  amplitude.  However,  the 
largest  errors  in  phase  are  most  often  associated  with  tidal  constituents  that  have  small  amplitudes. 
Thus,  it  is  thought  that  these  constituents  have  only  a  minor  impact  on  the  overall  accuracy  of 
computed  elevations.  In  some  instances,  particularly  with  high-frequency,  nonlinear  constituents, 
large  phase  errors  do  lead  to  increased  values  of  D.  The  large  phase  errors  exhibited  by  nonlinear 
tidal  components  are  not  surprising,  since  high-frequency  (and,  therefore,  shorter  wavelength)  tides 
exhibit  larger  spatial  gradients  in  phase.  Except  at  Christchurch  and  Lowestoft,  the  greatest  source 
of  error  as  measured  by  D  is  in  amplitude  of  the  semidiurnal  tidal  constituents.  These  constituents, 
as  noted  in  Sec.  5.0,  have  several  amphidromic  points  in  the  domain  and  it  is  well  known  that 
positioning  of  amphidromes  has  a  pronounced  effect  on  computed  phases. 

Consistently,  agreement  between  model-computed  elevations  and  the  data  is  far  better  at  those 
gauges  located  closest  to  the  open  boundaries.  This  outcome  is  not  surprising,  since  specified 
elevations  comprise  the  forcing  at  these  boundaries.  As  mentioned  previously,  the  quality  of  model 
predictions  at  Christchurch  may  be  strongly  influenced  by  its  proximity  to  an  amphidrome  and 
perhaps  by  omission  of  the  Isle  of  Wight.  No  correlation  between  error  and  gauge  location  are 
observed  with  respect  to  depths  at  the  gauge  location,  offshore  slope  at  the  gauge,  the  geographic 
location  (Continental  or  British  waters),  or  strength  of  the  tidal  nonlinearities  at  the  gauge. 


5.2  Velocity  Time  Series  Comparisons 

RMS  errors  from  comparisons  between  the  velocity  time  series  computed  by  ADCIRC-2DDI 
and  measured  velocity  time  series  at  each  of  the  eight  velocity  stations  is  presented  in  Table  3.  (The 
directional  errors  are  not  extremely  meaningful  as  the  largest  error  tends  to  occur  when  the  speed 
is  small  and  direction  is  changing  180°,  i.e.,  between  the  tidal  flood  and  ebb.)  Appendix  D  contains 
the  complete  set  of  time  series  plots  showing  these  comparisons.  From  the  plotted  comparisons 
shown  in  App.  D  and  from  velocity  comparisons  made  by  other  investigators  using  depth-integrated 
models,  it  seems  that  speed  observations  at  several  locations  are  of  marginal  value.  For  two  of  the 
gauges,  poor  correspondence  to  modeled  velocity  is  explained  by  the  unfortunate  vertical  placement 
of  the  gauges.  Gauges  4  and  8  are  located  near  the  bottom  at  a  depth  approximately  4/5  of  the  total 
water  depth  as  noted  by  Walters  (1987).  In  this  situation,  the  measured  speeds  are  expected  to  be 
much  lower  than  a  depth-averaged  speed  (e.g.,  Fig.  8a  and  b).  At  gauges  1,  3,  and  6,  comparisons 
to  ADCIRC-2DDI  computed  velocities  are  good  (e.g.,  Fig.  8c-e),  indicating  that  perhaps  at  these 
locations,  the  measured  velocities  are  well  represented  by  depth-averaged  velocities.  It  is  possible 
that  comparisons  to  current  measurements  of  this  type  could  be  improved  by  using  velocities 
computed  by  a  3D  tidal  model;  Lynch  and  Werner  (1984)  conducted  such  a  study  with  the  JONSDAP’76 
data  used  here,  but  no  improvement  was  observed. 
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Table  3  —  RMS  Errors  for  Raw  Velocity 
Model-Data  Comparisons 


Blain  and  Rogers 


GAUGE 

NUMBER 

SPEED 

RMS  ERROR  (m/s) 

DIRECTION 

RMS  ERROR  (deg) 

1 

.16 

12.7 

2 

.31 

19.9 

3 

.20 

39.9 

4 

.27 

51.9 

5 

.20 

26.8 

6 

.10 

40.2 

7 

.20 

37,7 

8 

.22 

31.4 

MEAN 

.21 

Encouragingly,  those  simulations  that  were  most  accurate  with  regard  to  elevation  also  compared 
best  with  respect  to  velocity.  It  should  be  noted  that  because  velocity  observations  are,  indeed,  the 
raw  data,  as  opposed  to  a  filtered  form  (i.e.,  the  elevation  observations),  there  is  no  reason  to 
filter  raw  model  velocity  computations  for  comparison. 


5.3  Performance  Comparison  to  Other  Models 

As  part  of  and  subsequent  to  the  Tidal  Flow  Forums,  a  sizable  number  of  researchers  have  used 
the  TFF  benchmark  case  to  assess  their  own  models’  performance  (e.g.,  Ozer  and  Jamart  1988; 
Jamart  and  Ozer  1989;  Yu  et  al.  1989;  Praagman  et  al.  1989a,  b;  Werner  and  Lynch  1987;  Gray 
et  al.  1987;  Walters  1987;  Lynch  and  Werner  1988;  Werner  and  Lynch  1989;  Gray  1989; 
Baptista  et  al.  1989;  Walters  and  Werner  1989;  Werner  1995).  RMS  errors  computed  for  the  last 
25  h  of  the  ADCIRC-2DDI  simulations  performed  here  are  compared  to  RMS  error  measures 
published  for  other  models.  In  general,  ADCIRC-2DDI  performed  well,  equaling  or  exceeding  the 
accuracy  of  other  models  for  both  raw  model  computations  and  filtered  prediction  comparisons 
based  on  the  11  TFF  constituents.  Table  4  contains  the  RMS  values  for  seven  models  including 
the  ADCIRC-2DDI  model  results  described  herein. 

These  inter-model  comparisons  should  be  approached  with  caution.  A  25-h  time  period  is 
considerably  short  for  meaningful  comparisons.  Comparisons  using  longer  time  series  may  lead  to 
an  entirely  different  ranking  of  model  “skill  level”  (Jamart  and  Ozer  1989).  Also,  variations  exist 
with  respect  to  the  computational  mesh  used  by  various  TFF  investigators.  Many  of  the  earlier 
modelers  used  the  original  FE  grid  discussed  previously.  Others  used  finite  difference  grids  and 
curvilinear  grids  that  offer  slightly  different  discretizations  of  the  domain.  Furthermore,  not  all 
modelers  compare  observations  to  the  model  solution  at  actual  gauge  locations;  rather,  several 
simply  compare  the  model  response  at  a  node  in  the  mesh  nearest  to  the  observation  location.  A 
minority  of  TFF  investigators  present  model  results  that  include  only  the  11  TFF  constituents; 
most  present  raw  model  time  series  (or  a  full  spectrum  solution  in  the  case  of  frequency  domain 
models). 


SPEED  (m/s)  SPEED  (m/s)  SPEED  (m/s) 
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DAY  (March  1976) 


O  OBSERVATIONS 
- BASE  SIMULATION  OUTPUT 


Fig.  8  —  Base  simulation  speeds;  velocities  at  gauges  (a)  4,  (b)  8,  and  (c)  1 
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Fig.  8  —  (cont.)  (d)  3  and  (e)  6 


Table  4  —  Model-Model  Comparisons  for  the  TFF  Benchmark  Case 


RMS  ERROR  (cm) 

RMS  ERROR  (cm) 

MODELS 

(FILTERED) 

(UNFILTERED) 

INVESTIGATORS 

ADCIRC-2DDI 

18.0 

32.5 

Blain  and  Rogers  1998 

WEQN 

22.0 

36.0 

Werner  and  Lynch  1989 

FADI 

42.0 

Yu  et  al.  1989 

WAQUA 

38.0 

Praagman  et  al.  1989a 

MU-model 

35.0 

Ozer  and  Jamart  1988 

MU-model 

27.0  (96  constituents) 

Jamart  and  Ozer  1989 

Walters 

28.0 

Walters  1987 
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6.0  SENSITIVITY  ANALYSES 
6.1  Description 

Sensitivity  analyses  are  conducted  to  investigate  tidal  dynamics  in  the  coastal  environment 
in  the  context  of  evaluating  the  implementation  and  performance  of  the  ADCIRC-2DDI  model  in 
simulating  these  flows.  Modularity  of  ADCIRC-2DDI  is  a  tremendous  asset  in  such  a  study  by 
allowing  terms  in  the  model  equations,  such  as  nonlinearities  and  various  forcing  mechanisms  to 
be  systematically  rendered  active  or  inactive. 

Four  sources  of  nonlinearity  are  contained  in  the  model  equations;  they  are  associated  with  the 
(1)  quadratic  bottom  stress  formulation,  which  is  dependent  on  V2,  (2)  finite  amplitude  terms,  those 
terms  involving  gradients  of  the  sea  surface  elevation,  V£,  (3)  spatial  derivatives  of  the  convective 
acceleration  terms,  e.g.,  v  •  V  v,  and  (4)  temporal  derivatives  of  the  convective  acceleration  terms, 
e.g.,  dV/dt.  “Switches”  contained  in  the  model  input  specifications  control  the  inclusion  or  exclu¬ 
sion  of  these  four  types  of  nonlinearity.  The  independent  effect  of  each  of  these  nonlinearities  is 
investigated  with  respect  to  their  influence  on  the  computed  tidal  dynamics  and  model  stability. 

Though  several  forms  of  forcing  at  the  open-ocean  boundary  can  be  implemented  in  ADCIRC, 
only  the  elevation  specified  condition  is  investigated  since  this  is  the  form  of  the  data  available  in 
the  TFF  benchmark.  Influence  of  the  frequency  content  of  this  boundary  elevation  forcing  on  the 
model-computed  response  is  examined.  Other  forcing  mechanisms  within  the  ADCIRC  model 
include  the  tidal  potential  forcing  that  serves  as  a  body  force  on  the  mass  of  the  water  body  and 
a  Coriolis  effect  accounting  for  latitude  of  the  modeled  region.  Additionally,  the  ADCIRC  model 
invokes  physical  and  numerical  parameters  such  as  the  bottom  friction  coefficient,  horizontal  mining 
coefficient,  model  timestep,  minimum  bathymetric  depth,  and  a  value  that  characterizes  the  form 
of  the  continuity  equation  invoked  x0.  All  forcing  and  parameters  just  mentioned  are  the  subject  of 
this  sensitivity  study  of  the  ADCIRC-2DDI  model. 


6.2  Nonlinearities 

6.2.1  Nonlinear  Bottom  Friction  Parameterization 

The  tidal  response  of  a  basin  is  strongly  influenced  by  topography.  Consequently,  the  bottom 
friction  parameterization  plays  a  significant  role  in  the  computed  tidal  response.  Bottom  friction 
within  ADCIRC-2DDI  can  be  represented  using  either  a  linear  relation  or  the  standard  quadratic 
parameterization  identified  in  the  presentation  of  Eqs.  (2.7)-(2.9).  For  the  nonlinear  bottom  friction 
formulation,  the  model  friction  factor  Cy  can  be  equated  to  physical  frictional  coefficients  such  as 
the  Manning  or  Chezy  forms.  In  contrast,  the  linear  bottom  friction  parameterization  has  a  friction 
factor  that  is  selected  based  on  considerations  of  mass  conservation,  model  stability,  and/or  model 
agreement  with  data  (i.e.,  calibration)  and  has  no  physical  basis. 

In  the  North  Sea/English  Channel,  both  linear  and  nonlinear  representations  of  the  bottom 
stress  clearly  capture  a  majority  of  the  tidal  signal  as  seen  at  four  stations  in  Fig.  9a-d.  No 
calibration  of  the  nonlinear  bottom  friction  factor  is  undertaken  for  these  simulations  so  that  con¬ 
sistency  with  the  TFF  benchmark  is  maintained.  A  nonlinear  bottom  friction  representation  is 
favored  because  of  its  physical  basis  for  selection  of  Cy.  Furthermore,  model  computations  compare 
slightly  better  with  observations  when  using  the  nonlinear  bottom  friction  representation  as  compared 
to  a  calibrated  linear  bottom  friction  formulation. 
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- NONLINEAR  BOTTOM 

FRICTION  ONLY 

- FULLY  LINEAR 

SIMULATION  (calibrated) 


Fig.  9  —  Effect  of  bottom  friction  parameterization;  tidal  gauges  (a)  Hoek  van  Holland  and  (b)  Walton 


A  linear  frictional  coefficient  of  6.0E-5  is  found  to  give  the  best  agreement  with  observations 
over  the  English  Channel/North  Sea.  The  use  of  much  larger  values  for  Cy  leads  to  significant 
damping  of  the  tidal  signal.  Simulations  applying  the  larger  linear  frictional  coefficients  are,  however, 
less  problematic  from  a  standpoint  of  mass  conservation  (to  be  discussed  further  in  Sec.  6.4.1). 

The  role  of  the  nonlinear  bottom  stress  on  different  tidal  frequencies  is  also  examined.  Nonlinearities 
associated  with  frictional  forcing  at  the  seabed  influence  to  a  greater  extent  the  semidiurnal  ( H 2) 
and  sextodiumal  (. H g)  constituents.  These  effects  are  illustrated  in  plots  of  global  amplitude  differences 
for  the  M2  and  Mg  constituents  that  contrast  simulations  employing  nonlinear  and  linear  bottom 
friction  formulations  (Fig.  10a  and  b). 


6.2.2  Nonlinear  Convective  Acceleration  Terms  ( Temporal  Derivatives) 

Westerink  et  al.  (1992b)  suggest  that  for  the  sake  consistency  between  the  continuity  and 
momentum  equations  and  for  better  mass  conservation  properties,  nonlinear  convective  acceleration 
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DAY  (March  1976) 

Fig.  9  —  (cont.)  (c)  Christchurch  and  (d)  Lowestoft 


terms  containing  time  derivatives  be  included  if  either  nonlinear  finite  amplitude  terms  or  spatial 
derivatives  of  the  nonlinear  convective  acceleration  terms  are  included  in  the  model  equations. 
Therefore,  a  sensitivity  analyses  for  the  time  derivative  terms  is  not  performed  independently  and 
no  conclusions  are  drawn  regarding  their  individual  effect  on  model  performance. 


6.2.3  Nonlinear  Finite  Amplitude  Terms 

The  total  water  depth  is  approximated  either  as  the  mean  water  depth,  H  =  h,  or  as  H  =  h  +  £, 
where  is  the  deviation  from  the  mean  water  depth  due  to  tidal  effects  (in  this  case).  The  latter 
representation  for  H  introduces  additional  nonlinearity  into  the  model  momentum  equations 
through  terms  that  multiply  water  surface  gradients  by  the  total  water  depth,  i.e.,  the  finite 

amplitude  terms.  Model  simulations  that  include  finite  amplitude  effects  in  addition  to  nonlinear 
bottom  stress  lead  to  slight  increases  in  predicted  sea  surface  height,  which  is  most  notable  at  the 
crests  and  troughs  of  the  tidal  cycle  (Fig.  lla-c). 
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15.0  15.5  16.0  16.5  17.0  17.5  18.0 

DAY  (March  1976) 

Fig.  11  —  Effect  of  nonlinear  finite  amplitude  terms;  tidal  gauges  at  (a)  Dover,  (b)  Boulogne,  and 

(c)  Calais 
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While  the  bottom  stress  parametrization  strongly  effects  the  semidiurnal  ( H2 )  and  sextodiurnal 
(H6)  tidal  frequencies,  nonlinear  finite  amplitude  terms  primarily  influence  the  quarter-diurnal  (//4) 
and  octodiurnal  (Hg)  constituents.  Global  amplitude  differences  computed  between  the  Af4  and  Mg 
tides  are  shown  in  Fig.  12a  and  b,  respectively,  which  contrast  simulations  with  and  without  the 
finite  amplitude  effects. 


Fig.  12  —  Global  amplitude  differences  (cm);  amplitude  with  finite  amplitude 
terms  -  amplitude  without  finite  amplitude  terms,  (a)  M4  and  (b)  M8 
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Differences  in  the  tidal  constituent  between  the  two  simulations  are  considerable  west  of 
the  Cherbourg  peninsula  and  near  the  Dutch  coastline.  A  time  series  comparing  the  quarter-diurnal 
constituents  from  the  two  simulations  to  observations  at  a  gauge  location  on  the  Dutch  coast  is 
shown  in  Fig.  13. 

6.2.4  Nonlinear  Convective  Acceleration  Terms  ( Spatial  Derivatives) 

Inclusion  of  the  spatial  derivatives  of  the  convective  acceleration  terms  in  the  model  equations 
together  with  nonlinear  bottom  stress  terms  results  causes  a  slight  reduction  in  the  surface  water 
elevation  at  the  crests  and  troughs  of  the  tidal  cycles  in  Fig.  14a  and  b.  A  notable  difference  in 
the  computed  tides  is  observed  when  both  the  finite  amplitude  and  convective  acceleration  terms 
are  either  included  or  excluded  from  the  model  equations  (see  Fig.  15a— d  comparing  the  full 
nonlinear  solution  to  a  solution  that  has  bottom  stress  as  the  only  nonlinearity). 

The  influence  of  the  nonlinear  mechanisms  associated  with  the  finite  amplitude  and  convective 
terms  is  clearly  evident  in  the  contributions  of  five  nonlinear  TFF  constituents  to  the  total  tidal 
elevation  as  shown  in  Fig.  16a— c.  These  five  constituents  are  largest  when  all  nonlinear  terms  are 
included  in  the  model  computations.  Furthermore,  the  strength  of  these  nonlinear  tidal  constituents 
indicates  that  nonlinear  interactions  within  the  domain  (facilitated  through  the  nonlinear  terms  in 
the  model  equations)  as  opposed  to  nonlinear  forcing  at  the  boundary  are  largely  responsible  for 
generation  of  the  nonlinear  tidal  response.  This  is  most  evident  at  Boulogne,  which  is  far  from 
open  boundaries  where  the  five  nonlinear  constituents  are  forced.  The  North  Sea/English  Channel 
is  an  excellent  example  of  nonlinear  tidal  generation  in  shallow  waters  and  demonstrates  the 
necessity  of  a  nonlinear  hydrodynamic  model  for  the  prediction  of  shallow-water  tidal  dynamics. 

Not  surprisingly,  tidal  circulation  is  also  altered  by  the  inclusion  or  omission  of  the  finite 
amplitude  and  convective  acceleration  terms.  The  speed  of  the  flow  is  particularly  affected  with 
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Fig.  13  —  Effect  of  finite  amplitude  terms  on  three  quarter-diurnal  TFF  constituents;  tidal  gauge 

at  Hoek  van  Holland 
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DAY  (March  1976) 

Fig.  14  —  Effect  of  convective  acceleration  terms,  tidal  gauges  at  (a)  Cherbourg  and  (b)  Lowestoft 


directions  altered  to  a  much  lesser  degree.  Global  differences  computed  with  respect  to  speed 
emphasize  this  point  and  are  shown  in  Fig.  17  for  simulations  with  and  without  the  finite  amplitude 
and  convective  acceleration  terms  collectively. 

6.2.5  Effect  of  Nonlinearities  on  Accuracy  of  Filtered  Model  Solutions 

Here,  the  inclusion  of  model  nonlinearities  on  filtered  model-to-data  comparisons  is  examined. 
The  model  predictions  used  for  comparison  contain  only  those  11  TFF  constituents  included  in 
the  observations.  Table  5  contains  proportional  errors  computed  from  comparisons  between  eleva¬ 
tion  time  series  for  a  3-d  time  period  (0000  15-18  Mar).  Predictions  from  10  of  the  11  tidal 
elevation  gauges  are  included  in  computations  of  the  average  error  in  each  time  series;  Christchurch 
is  excluded  for  reasons  previously  discussed.  Clearly,  when  nonlinear  terms  are  active  in  the  model 
equations,  model-to-data  comparisons  significantly  improve  as  reflected  by  smaller  proportional 


Coastal  Tide  Prediction 


35 


Fig.  15  • 


15.5  16.0  16.5  17.0  17.5  18.0 

DAY  (March  1976) 

•  Effect  of  finite  amplitude  and  convective  acceleration  terms,  tidal  gauges  at  (a)  Lowestoft 

and  (b)  Calais 


errors.  Figure  18a-d  further  illustrates  the  effect  of  full  nonlinearity  on  the  computed  tidal  response 
as  compared  to  the  measured  sea  surface  heights  at  four  gauges.  Recall  that  Fig.  9a-d  compared 
the  model  response  from  the  nonlinear  bottom  friction  alone  to  the  field  observations. 


6.3  Forcing 

6.3.1  Inclusion  of  Tidal  Constituents 

For  this  study  case,  six  linear  and  five  nonlinear  tidal  constituents  are  provided  at  the  open 
boundaries.  Often,  tidal  dynamics  on  the  continental  shelf  are  difficult  to  characterize.  In  particular, 
nonlinear  tidal  constituents  are  generally  not  known  a  priori  so  that  most  often,  only  deep-water, 
primary  tidal  constituents  are  applied  as  boundary  forcing.  For  the  North  Sea,  two  forms  of  the 
boundary  forcing  are  examined:  the  first  contains  only  six  primary  TFF  constituents  and  the  second 
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O  OBSERVATIONS 
-  -  FULLY  NONLINEAR  MODEL 

- NONLINEAR  BOTTOM 

FRICTION  ONLY 


Fig.  15  —  (cont.)  (c)  Hoek  van  Holland  and  (d)  Dover 


uses  all  11  TFF  constituents.  Inclusion  of  the  nonlinear  tidal  constituents  does  have  a  slight  effect 
on  model  results  as  seen  in  Fig.  19a  and  b.  The  differences  between  tidal  predictions  forced  with 
and  without  nonlinear  tides  is  at  least  partially  a  result  of  the  open-ocean  boundary  location  in 
relatively  shallow  water  (<120  m)  where  nonlinear  tides  are  significant.  One  would  expect  errors 
for  all  tidal  constituents  if  any  other  than  the  11  TFF  constituents  are  significant  at  the  open 
boundaries,  due  to  nonlinear  interactions  with  improperly  forced  non-TFF  tidal  constituents  (Werner 
and  Lynch  1989). 

Tables  6a  and  b  list  the  calculated  values  for  the  error  in  the  complex  plane  D  (averaged  for 
11  gauges)  for  each  of  the  11  TFF  constituents  computed  by  the  model.  Comparison  of  the  simulations, 
which  include  all  constituents  in  the  boundary  forcing  (Table  6a)  and  those  using  only  six  primary 
tidal  constituents  for  forcing  (Table  6b),  clearly  demonstrates  the  benefit  of  including  nonlinear 
constituents  in  the  boundary  forcing. 
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O  OBSERVATIONS 

-  NONLINEAR  BOTTOM 

FRICTION  ONLY 
- FULLY  NONLINEAR 


Fig.  16  —  Effect  of  finite  amplitude  and  convective  acceleration  terms  on  five  nonlinear  TFF 
constituents,  tidal  gauges  at  (a)  St.  Malo,  (b)  Boulogne,  and  (c)  Hoek  van  Holland 
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Fig.  17  —  Global  speed  differences  (m/s);  speed  with  finite  amplitude  and  convective  acceleration 
terms  -  speed  without  finite  amplitude  and  convective  acceleration  terms 


Table  5  —  The  Role  of  Nonlinearity  on  Model-Data  Comparisons 


NONLINEAR 

NONLINEAR 

NONLINEAR 

CONVECTIVE 

CONVECTIVE 

DIMENSIONLESS 

FINITE 

ACCELERATION 

ACCELERATION 

PROPORTIONAL 

NONLINEAR 

AMPLITUDE 

(Spatial 

(Temporal 

ERROR, 

BOTTOM  STRESS 

TERMS 

Derivatives) 

Derivatives) 

£2  xlOO 

yes 

yes 

yes 

yes 

1.07 

yes 

yes 

no 

yes 

1.16 

yes 

no 

yes 

yes 

1.82 

yes 

no 

no 

no 

2.93 

no 

no 

no 

no 

3.99 

(with  calibration) 
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15.0  15.5  16.0  16.5  17.0  17.5  18.0 

DAY  (March  1976) 

Fig.  18  —  Effect  of  finite  amplitude  and  convective  acceleration  terms  on  filtered  model  output,  tidal 
gauges  at  (a)  Boulogne  and  (b)  Hoek  van  Holland 


6.3.2  Tidal  Potential  and  Variable  Coriolis  Forcing 

To  include  tidal  potential  forcing  (a  body  force  on  the  water  mass  due  to  the  gravitational 
forces  of  the  sun  and  the  moon)  and  variable  Coriolis  forcing,  a  transformation  of  the  model  domain 
is  made  from  Cartesian  coordinates  (supplied  by  the  TFF  organizers)  to  spherical  coordinates: 


X  =  + 


(6.1) 


<t>  =  bo  7 .  (6.2) 

Here,  x  and  y  are  Cartesian  coordinates  referenced  to  a  longitude  and  latitude  origin  near  the  center 
of  the  grid  (kQ,  and  r  is  the  radius  of  the  Earth.  For  purposes  of  computing  a  variable  Coriolis 
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O  OBSERVATIONS 
--- . FULLY  NONLINEAR  MODEL 

- NONLINEAR  BOTTOM 

FRICTION  ONLY 


Fig.  18  —  (cont.)  (c)  Dover  and  (d)  Christchurch 


forcing,  the  centroid  of  the  domain  is  selected  at  0.3°  latitude  and  50.75°  longitude.  In  simulations 
including  either  the  tidal  potential  or  the  variable  Coriolis  forcing,  the  added  complexity  does  not 
result  in  any  significant  change  in  model  tidal  response.  At  some  locations,  velocity  directions  are 
slightly  varied  with  the  inclusion  of  tidal  potential  forcing  and  variable  Coriolis  forcing  (c-g-, 
Fig.  20a-c),  but  the  effects  are  inconclusive.  Though  tidal  potential  forcing  and  the  variability  of 
Coriolis  forcing  seem  to  be  insignificant  for  this  relatively  small  study  area,  a  larger  domain 
covering  a  wider  range  of  latitude  will  necessitate  the  inclusion  of  these  two  forcing  mechanisms. 


6.4  Parameters 

6.4.1  Weighting  Factor 

The  GWCE  weighting  factor,  x0,  found  in  Eqs.  (2.19)  and  (2.20),  allows  flexibility  in  the 
formulation  of  the  equations  implemented  by  ADCIRC-2DDI.  For  x0  =  0,  the  GWCE  becomes  a 
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15.0  15.5  16.0  16.5  17.0  17.5  18.0 

DAY  (March  1976) 


O  OBSERVATIONS 

- BOUNDARY  FORCING  BY 

6  LINEAR  CONSTITUENTS 

- BOUNDARY  FORCING  BY 

11  CONSTITUENTS 


Fig.  19  —  Effect  of  constituents  in  boundary  forcing,  tidal  gauges  at  (a)  Hoek  van  Holland 

and  (b)  Dover 


pure  wave  equation,  whereas  large  values  of  xa  shift  the  GWCE  toward  a  primitive  form  of  the 
continuity  equation.  Westerink  et  al.  (1992b)  state  that  xQ  must  be  carefully  chosen,  as  overly  small 
values  lead  to  mass  conservation  problems  and  too-large  values  result  in  instabilities  caused  by 
2  •  Ax  waves.  It  is  well  known  that  FE  approximations  of  the  primitive  equations  suffer  from 
spurious  oscillations;  this  was  a  common  problem  for  early  generation  FE  models. 

For  cases  in  which  nonlinear  bottom  stress  is  included,  Westerink  et  al.  (1992b)  suggest  setting 
the  weighting  factor  equal  to  the  maximum  x*  for  the  simulation  (e.g.,  x*  =  Cf(u2  +  v2)max/hmin) 
and  the  model  friction  factor  equal  to  the  bottom  friction  coefficient  Cf.  Our  experience  indicates 
that  x*  =  [Cf(u2  +  v2)/h]max  is  less  conservative,  but  nonetheless  acceptable.  When  a  linear  bottom 
stress  formulation  is  applied,  it  is  suggested  that  x0  and  the  model  friction  factor  be  assigned 
identical  values.  Typically,  xQ  is  determined  by  calibration  or  is  based  on  the  user’s  judgment  and 
experience  for  a  given  application.  Mass  conservation  is  intimately  tied  to  the  specification  of  the 
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Table  6a  —  Elevation  Errors  for  Boundary  Forcing  with  11  TFF  Constituents 


GAUGE 

LOCATION 

01 

K1 

M2 

S2 

N2 

K2 

M4 

MS4 

MN4 

M6 

2MS6 

MEAN 

St.  Malo 

.017 

.018 

.180 

m 

m 

m 

m 

.069 

.044 

.018 

.018 

.057 

Cherbourg 

.015 

.016 

.107 

E 

1 

m 

in 

.017 

.014 

.012 

.031 

Dieppe 

.015 

.027 

.170 

m 

ESI 

.048 

.035 

.037 

.025 

.013 

.021 

.053 

Boulogne 

.014 

.026 

.100 

.072 

.040 

.026 

.032 

.031 

.020 

.014 

.014 

.035 

Calais 

.015 

.015 

.242 

.118 

.068 

.025 

.054 

.059 

.020 

.035 

.024 

.061 

Zeebrugge 

.011 

.040 

.120 

.092 

.138 

.025 

.041 

.016 

.024 

.020 

.039 

.052 

Hoek  van  Holland 

.025 

.006 

.109 

.014 

.006 

.011 

.078 

.018 

.026 

.019 

.022 

.030 

Walton 

.025 

.035 

.223 

.105 

.103 

.027 

.020 

.029 

.023 

.017 

.012 

.056 

Dover 

.003 

.016 

.092 

.034 

.031 

.017 

.073 

.034 

.058 

.017 

.019 

.036 

Christchurch 

.013 

.004 

.045 

.033 

.036 

.008 

.078 

.041 

.023 

.071 

.060 

.037 

Lowestoft 

.016 

.043 

.040 

.018 

.025 

.010 

.023 

.015 

.002 

.025 

.061 

.025 

MEAN 

.015 

.022 

.013 

.066 

.056 

.023 

.052 

.033 

.025 

.023 

.027 

.043 

Table  6b  —  Elevation  Errors  for  Boundary  Forcing  with  Six  TFF  Constituents 


GAUGE 

LOCATION 

Ol 

K1 

M2 

S2 

N2 

K2 

M4 

MS4 

MN4 

M6 

2MS6 

MEAN 

St.  Malo 

.017 

.200 

.238 

.151 

■ m 

.038 

mm 

.088 

Cherbourg 

.015 

.103 

|5j 

.101 

.058 

El 

.024 

m 

.044 

Dieppe 

.016 

.163 

.114 

jlll 

.064 

.047 

RSI 

.008 

EH 

.055 

Boulogne 

.015 

.097 

.068 

.038 

Ell 

.090 

.072 

.041 

.012 

.029 

.047 

Calais 

.016 

.015 

.242 

.118 

.066 

.024 

.065 

.049 

.029 

.027 

.004 

.060 

Zeebrugge 

.012 

.041 

.120 

.093 

.136 

.025 

.046 

.007 

.021 

.011 

.037 

.050 

Hoek  van  Holland 

.026 

.007 

.106 

.012 

.004 

.011 

.115 

.057 

.028 

.052 

.050 

.043 

Walton 

.026 

.035 

.225 

.107 

.102 

.027 

.030 

.008 

.018 

.010 

.013 

.055 

Dover 

.003 

.016 

.095 

.033 

.029 

.017 

.122 

.080 

.069 

.016 

.031 

.046 

Christchurch 

.013 

.004 

.042 

.034 

.036 

.008 

.085 

.040 

.031 

.073 

.063 

.039 

Lowestoft 

.017 

ESI 

.040 

.017 

.024 

.010 

.031 

.028 

.016 

.019 

.016 

.024 

MEAN 

.016 

m 

.013 

.065 

.054 

.023 

.090 

.054 

.038 

.026 

.030 

.050 

parameter  rQ.  Mass  balance  properties  are  known  to  be  suspect  as  the  form  of  the  GWCE  continuity 
equation  shifts  towards  a  pure  wave  equation;  for  this  situation,  the  largest  errors  have  been  shown 
to  be  located  at  the  open  boundaries. 

Mass  balance  properties  of  model  simulations  are  investigated  here  to  gain  better  insight  into 
the  selection  of  t0  and  determine  the  influence  of  other  factors  on  mass  conservation.  The  mass 
accumulation  over  time  is  shown  in  Fig.  21a  and  b  for  simulations  using  a  linear  bottom  friction 
parameterization.  Since  the  volume  change  in  mass  and  mass  flux  into  the  domain  should  be  equal, 
mass  balance  errors  are  represented  as  the  difference  between  these  two  curves.  It  is  uncertain 
whether  the  smaller  discrepancy  in  mass  for  the  uncalibrated  linear  friction  simulation  is  a  result 
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DAY  (March  1976) 


O  OBSERVATIONS 

-  WITHOUT  TIDAL  POTENTIAL 

FORCING  AND  VARIABLE 
CORIOLIS  FORCING 

- WITH  TIDAL  POTENTIAL 

FORCING  AND  VARIABLE 
CORIOLIS  FORCING 


Fig.  20  —  Effect  of  variable  Coriolis  forcing,  velocity  gauges  (a)  2,  (b)  7,  and  (c)  8 


CUMULATIVE  MASS  (m3)  CUMULATIVE  MASS  (m3)  CUMULATIVE  MASS  (m3) 
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Fig.  21  —  Mass  balance  analysis,  linear  models,  (a)  calibrated  (tauO  =  ffactor  =  6E-5), 
(b)  uncalibrated  (tauO  =  6E-5;  ffactor  =  Cf  =  6E-4),  and  (c)  fully  nonlinear  model 
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of  the  small  difference  between  the  model  parameters  x0  and  FFACTOR,  the  linear  friction 
coefficient,  or  whether  it  is  merely  due  to  the  weaker  tidal  response  caused  by  the  higher  fric¬ 
tional  coefficient.  Mass  balance  is  significantly  improved  when  using  a  nonlinear  bottom  stress 
representation  (Fig.  21c)  as  compared  to  the  model  simulations  with  calibrated  linear  bottom  friction 
(Fig.  21a).  For  all  the  simulations,  water  volume  changes  are  found  to  be  reasonable  while  the  mass 
flux  is  consistently  in  error. 

Checks  of  mass  conservation  are  also  conducted  for  three  subsections  of  the  grid:  the  English 
Channel  constriction,  a  region  along  the  French  coast,  and  the  western  open  boundary  (see 
Fig.  22a-c).  For  the  French  coastal  section  and  the  English  Channel  constriction,  mass  conservation 
computations  mimic  those  produced  for  the  entire  grid.  The  mass  flux  into  the  domain  remains 
high,  skewing  mass  balance  for  the  computed  solution.  As  expected,  Fig.  22c  shows  large  mass 
balance  errors  occur  at  the  open  boundary.  Contrary  to  the  interior  locations  in  the  domain,  the 
mass  flux  out  of  this  region  is  greater  than  the  volume  of  mass  lost. 

The  mass  conservation  properties  of  the  GWCE  equation  continue  to  be  the  subject  of  current 
research,  e.g.,  Lynch  and  Holboke  (1997),  Westerink  et  al.  (1994b). 

6.4.2  Time-Step  Size 

Time-step  size  selection  is  based  only  on  considerations  of  stability  through  the  computation 
of  an  acceptable  Courant  number  (Lapidus  and  Pinder  1982).  For  example,  a  timestep  appropriate 
to  the  original  TFF  coarse  grid  is  determined  as: 


A t  <  300s  = 


1.5  Axmin 

\J  S^1  max 


1.5  (9  km) 

\/  9.81  (m/s2)  (110  m) 


(6.3) 


A  Courant  number  of  1.5  is  suggested  by  Westerink  et  al.  (1992b),  though  experience  suggests  that 
values  less  than  1.0  are  optimal.  Experiments  using  a  variety  of  timesteps  demonstrate  that  as  long 
as  the  stability  criterion  is  met,  time-step  size  has  no  observed  influence  on  model  computations. 


6.4.3  Lateral  Eddy  Viscosity 

A  sensitivity  analysis  on  the  lateral  mixing  term  is  performed  by  selecting  eddy  viscosity 
coefficient  values  in  the  range  from  0  to  100  m2/s,  three  orders  of  magnitude.  No  significant  effect 
of  the  lateral  eddy  viscosity  on  model-computed  tides  is  observed.  The  horizontal  mixing  terms 
can  be  expected  to  play  a  more  important  role  when  the  environment  is  advection  dominated  or 
small-scale,  nearshore  flows  are  considered. 

6.4.4  Minimum  Depth 

All  model  simulations  to  this  point  have  not  included  a  mechanism  for  shoreline  inundation. 
For  these  simulations,  it  is  determined  that  a  minimum  depth  of  approximately  8  m  is  required 
when  nonlinear  finite  amplitude  terms  are  included  in  the  model  formulation.  The  use  of  smaller 
minimum  depths  (e.g.,  6.0  m)  results  in  negative  depths  indicative  of  drying  of  the  grid  elements. 
Because  no  “standard”  minimum  depth  is  suggested  by  the  TFF  organizers,  a  value  of  Ha  =  10  m 
is  specified  in  all  non-wetting/drying  simulations,  including  those  simulations  that  do  not  involve 
nonlinear  finite  amplitude  terms  for  the  sake  of  consistency. 


46 


Blain  and  Rogers 


xIO5 


-3-2-10  1  2  3 

x(m)  xIO5 


—  DOMAIN  BOUNDARY 

o  NODES  IN  THE  ENGLISH 
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-  NET  FLUX 

-  ACCUMULATION 


Fig.  22  —  (a)  Subgrid  sections  used  in  mass  balance  analysis;  mass  balance  analysis  fully  nonlinear  models  and 

(b)  English  Channel  constriction 


For  those  model  simulations  that  employ  a  shoreline  inundation  mechanism,  the  specification 
of  a  minimum  depth  is  not  necessary. 


6.5  Implementation  of  Wetting  and  Drying 

Within  ADCIRC-2DDI,  shoreline  inundation  is  handled  through  a  numerical  wetting  and  drying 
procedure  in  which  mesh  elements  can  become  active  or  inactive  throughout  the  duration  of  the 
simulation.  This  wetting  and  drying  feature  in  ADCIRC-2DDI  is  implemented  successfully  for 
the  North  Sea/English  Channel  study  case.  Allowing  the  wetting  and  drying  of  elements  makes 
possible  the  inclusion  of  nonlinear  finite  amplitude  terms  without  imposing  a  large  minimum  depth. 
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TIME  (days  into  simulation) 

Fig.  22  —  (cont.)  (c)  French  coast  and  (d)  western  open  boundary 


During  model  simulations  of  the  North  Sea/English  Channel  when  the  wetting  and  drying  feature 
is  implemented,  drying  occurs  at  many  of  the  11  coastal  elevation  gauge  locations  preventing  any 
model-data  comparisons  based  on  elevation.  A  notable  observation  is  that  drying  occurs  at  coastal 
nodes  with  relatively  large  bathymetric  depths  (e.g.,  8  m  below  mean  sea  level)  with  an  indication 
that  this  drying  is  due  to  land-locking.  Based  on  the  location  of  such  a  node  in  Fig.  23,  land-locking 
seems  unrealistic. 

Mass  conservation  improves  when  wetting  and  drying  is  active  as  indicated  by  the  global  mass 
flux  and  mass  accumulation  curves  shown  in  Fig.  24a  and  b.  The  reason  for  this  improvement  is 
related  to  the  elimination  of  a  minimum  depth  constraint  that  serves  as  a  source  of  mass  to  the 
system.  Furthermore,  when  using  the  wetting/drying  feature,  alterations  made  to  the  numerical 
solution  procedure  may  contribute  toward  improved  mass  conservation. 
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Fig.  23  —  A  suspect  example  of  ADCIRC’s  wetting  and  drying  routine 
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24  —  Mass  balance  analysis,  fully  nonlinear  model,  wetting  and  drying,  (a)  without 

and  (b)  with 
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Wetting  and  drying  is  a  very  recent  feature  of  ADCIRC-2DDI,  and  by  the  admission  of  its 
developers,  is  largely  untested  with  needed  improvements  acknowledged.  The  option  of  element 
wetting  and  drying  does  offer  possibilities  for  additional  realism  and  further  model  stability,  which 
in  turn  permits  the  inclusion  of  all  model  nonlinearities  in  simulations  using  ADCIRC-2DDI. 

It  should  be  noted  that  implementation  of  the  wetting  and  drying  feature  does  not  require 
specification  of  a  special  land  boundary  type,  such  as  a  weir  boundary.  Specification  of  this  type 
of  land  boundary  is  only  necessary  when  representing  structures  such  as  alongshore  dikes  or  seawalls. 


7.0  MESH  RESOLUTION  ISSUES 

Mesh  resolution  resulting  from  the  placement  of  discrete  nodes  throughout  the  domain  cannot 
be  entirely  separated  from  the  model-computed  response.  Valuable  insight  into  an  appropriate  grid 
resolution  for  a  given  problem  is  gained  by  considering  the  scale  of  the  physical  processes  being 
modeled;  in  this  case,  the  physical  scales  are  associated  with  tidal-induced  circulation.  Generally 
and  specifically  in  the  English  Channel/North  Sea,  tidal  forcing  produces  very  long  waves  (e.g., 
semidiurnal  tides  are  O  (1000  km))  that  are  readily  resolved  by  a  relatively  coarse  mesh  having 
nodal  spacing  of  ~10  km.  However,  irregular  land  boundaries  whose  variability  occurs  at  smaller 
spatial  scales,  O  (less  than  10  km),  induce  large  velocity  gradients.  Proper  computation  of  the  tidal 
circulation  in  these  circumstances  requires  significantly  more  mesh  refinement  than  is  needed  to 
predict  sea  surface  height. 

Mesh  resolution  also  directly  influences  truncation  errors  associated  with  the  discrete  equations 
that  can  lead  to  degradation  of  the  computed  solution  (i.e.,  Westerink  et  al.  1994c;  Hagen  and 
Westerink  1995,  1996;  Westerink  and  Roache  1996).  The  effects  of  resolution  on  both  the  computed 
physics  and  the  numerical  approximation  is  determined  through  a  grid  convergence  study  (e.g., 
Dietrich  et  al.  1990;  Lardner  and  Song  1992;  Blain  et  al.  1998).  As  a  minimum,  grid  convergence 
studies  compare  solutions  computed  over  two  meshes  of  different  resolution.  If  the  modeled  physics 
are  not  adequately  resolved  over  a  coarse  grid,  the  computed  coarse  grid  solution  will  differ 
significantly  from  the  fine  grid  result.  An  infinite  degree  of  resolution  will  remove  inadequacies 
due  to  mesh  refinement,  but  this  approach  is  clearly  not  practical.  A  practical  objective  is  to  select 
a  resolution  that  captures  the  physics  adequately  and  is  computationally  efficient.  Figure  25  illustrates 
this  motivation  for  a  grid  convergence  study. 


7.1  Computation  of  Grid  Errors 

Model  simulations  over  the  North  Sea/English  Channel  are  conducted  for  three  discrete  meshes 
with  varying  levels  of  nodal  point  resolution.  The  three  grids  are  referred  to  as:  G9  with  nodal 
spacing  ranging  from  9-15  km,  G4  with  nodal  spacing  ranging  from  4-8  km,  and  G2  with  nodal  spacing 
ranging  from  2-4  km.  Properties  of  these  grids  are  described  further  in  Table  7.  Also  included  is 
an  indication  of  the  typical  computation  time  for  a  45-d  simulation  on  a  Sun  UltraSparc  workstation. 

Figure  26a  and  b  contains  times  series  plots  at  two  locations  near  the  coast  of  the  Cherbourg 
peninsula  shown  in  Fig.  27  for  simulations  using  the  grids  G4  and  G9;  note  that  grid  G4  is  simply 
a  double  refinement  of  grid  G9.  At  two  locations,  node  192  and  node  1605,  elevation  differences 
are  typically  large  when  checked  against  other  locations  in  the  domain.  Overprediction  at  one  node 
together  with  the  corresponding  underprediction  at  a  nearby  node  is  thought  to  be  due  to  increased 
velocity  gradients  in  that  area.  Figure  28  depicts  notable  changes  in  the  velocity  magnitude  computed 
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ERROR  DUE  TO  INADEQUATE  RESOLUTION 
Fig.  25  —  Optimization  of  grid  resolution 


Table  7  —  Properties  of  Computational  Grids  for  the  North  Sea/English  Channel 


GRID 

NODAL  SPACING 
(km) 

NUMBER  OF 
NODES 

NUMBER  OF 
ELEMENTS 

CPU  FOR  45 -d 
SIMULATION  (h) 

G9 

9-15 

911 

1613 

0.7 

G4 

4-8 

3434 

6452 

3.5 

G2 

2-4 

13319 

25808 

35.3 
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Fig.  26  —  Elevation  time  series  near  Cherbourg,  (a)  node  192  and  (b)  node  1605 


Fig.  27  —  Cherbourg  Peninsula,  shown  with  FE  grid 
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DAY  (March  1976) 


—  GRID  G9 
GRID  G4 


Fig.  28  —  Speed  comparisons  at  node  192 


over  grids  G9  and  G4  at  node  192.  Differences  in  the  computed  circulation  patterns  are  likely  due 
to  changes  in  spatial  bathymetry  gradients  that  are  a  direct  consequence  of  the  differing  mesh 
resolutions. 

Elevation  time  series  at  node  192  are  reconstructed  using  several  combinations  of  tidal  constituents 
as  computed  over  grids  G9  and  G4.  Comparisons  of  these  time  series  are  shown  in  Fig.  29a-e  for 
fully  nonlinear  simulations.  Note  that  the  coarse  grid  results  have  been  interpolated  onto  the  fine 
grid  to  facilitate  difference  computations.  Tidal  amplitudes  based  on  either  all  11  TFF  constituents 
(Fig.  29a)  or  the  six  astronomical  TFF  constituents  (Fig.  29b)  are  larger  when  computed  over  the 
coarser  mesh.  Figure  27c  shows  notable  differences  in  computation  of  the  five  TFF  nonlinear  tides 
over  the  two  meshes.  The  inclusion  of  56  tidal  constituents  (which  approximates  the  raw  model 
solution)  in  the  time  series  reconstruction  (Fig.  29d)  results  in  deeper  troughs  of  the  tidal  signal 
over  the  coarse  grid.  Dynamics  at  the  trough  of  a  tidal  wave  are  such  that  current  speeds  increase 
due  to  a  reduced  water  depth  at  the  trough.  Higher  velocity  magnitudes  and  sharp  velocity  gradients 
near  the  trough  may  not  be  well  represented  over  the  coarse  mesh  and  could,  in  fact,  lead  to 
degradation  of  the  velocity  solution.  For  the  tidal  signal  reconstruction  based  on  the  45  non-TFF 
constituents  (Fig.  29e),  elevations  are  consistently  less,  using  the  coarse  grid  solution  as  compared 
to  the  fine  grid  solution.  Differences  between  model  solutions  over  each  grid  are  thought  to  be 
primarily  caused  by  alterations  in  the  circulation  pattern  due  to  bathymetry  deviations  at  different 
resolutions. 

Domain-wide  tidal  elevations  computed  over  the  coarse  grid,  G9,  are  now  compared  to 
domain-wide  elevations  computed  over  the  fine  mesh  grid,  G4,  producing  difference  errors  between 
the  two  mesh  solutions.  Figure  30a  and  b  present  histograms  of  elevation  differences  between  each 
pair  of  grid  simulations,  i.e.,  G4  and  G9,  G2  and  G4,  and  represent  a  composite  for  the  entire 
domain.  The  average  difference  between  the  grid  G9  and  G4  solutions  is  relatively  small,  approxi¬ 
mately  1.4  cm,  and  errors  are  spread  in  a  Gaussian-like  distribution  over  errors  ranging  to  ±0.05  m. 
In  contrast,  elevations  computed  using  grids  G4  and  G2  exhibit  much  more  similarity  than  compari¬ 
sons  between  solutions  over  grids  G9  to  G4.  In  Fig.  30b,  the  average  magnitude  of  elevation 
differences  between  solutions  computed  over  grids  G4  and  G2  is  0.4  cm,  a  considerable  reduction 
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DAY  (March  1976) 


Fig.  29  —  Tidal  elevations  at  node  192,  (a)  11  TFF  constituents  only,  (b)  six  astronomical  TFF 
constituents  only,  and  (c)  five  nonlinear  TFF  constituents  only 
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Fig.  29  —  (cont.)  (d)  56  constituents  and  (e)  45  non-TFF  constituents  only 


ELEVATION  DIFFERENCE  (m) 


Fig.  30  —  Grid  error  histograms,  simulations  with 
nonlinear  bottom  friction  and  nonlinear  finite 
amplitude  terms,  (a)  G4  and  G9  and  (b)  G2  and  G4 
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over  the  1.4  cm  average  for  the  G9  to  G4  comparison.  Furthermore,  nearly  the  entire  distribution 
of  error  is  within  ±0.01  m.  Resolution  effects  are  quickly  reduced  as  refinement  increases,  indicating 
that  the  computed  model  solution  is  rapidly  converging  toward  the  asymptotic  solution,  which 
would  result  for  infinite  grid  resolution.  Westerink  and  Roache  (1996)  suggest  a  test  to  determine 
whether  a  grid  falls  within  the  asymptotic  range  using  a  Richardson-based  error  estimator  for  the 
fine  and  coarse  grid  solutions. 

A  Richardson-based  error  measure  derives  an  error  estimate  for  model  solutions  computed  over 
both  a  coarse  and  fine  mesh  (see  Blain  et  al.  1998  for  implementation).  For  the  coarse  mesh,  the 
Richardson-based  error  estimate  is  cast  as  Ec  and  for  the  fine  mesh,  the  error  is  represented  as  EF. 
Expressions  for  these  errors  are  derived  by  Roache  (1994)  and  presented  below  for  this  application 
in  which  the  fine  grid  is  obtained  by  doubling  the  mesh  refinement  of  the  coarse  grid  and  the  spatial 
discretization  of  the  discrete  model  equations  is  second  order: 


(6.4) 


EF=\  W-  (6.5) 

Note  that  £max  is  the  maximum  elevation  difference  between  the  coarse  and  fine  grid  solutions  at 
a  specified  instant  in  time.  Normalized  error  estimates  Ec/n,  EF/n  are  obtained  by  dividing  the  error 
estimates  of  Eqs.  (6.4)  and  (6.5)  by  the  maximum  elevation  over  the  domain  during  the  time 
interval  of  interest. 

To  determine  if  a  grid  lies  in  the  asymptotic  range,  a  simple  comparison  between  fine  and 
coarse  grid  error  estimates,  Ec,  EF  for  that  grid  is  made.  If  Ec  and  EF  for  that  grid  yield  similar 
magnitudes,  then  the  grid  is  assumed  to  lie  in  the  asymptotic  range.  Figure  31a  and  b  shows  results 
from  these  comparisons  for  two  simulations,  both  using  a  nonlinear  bottom  friction  parameterization 
and  one  including  the  finite  amplitude  terms.  Slight  differences  between  Ec  and  EF  indicate  that 
G4  does  not  lie  in  the  asymptotic  range.  However,  these  differences  are  not  great,  suggesting 
that  G4  is  near  the  asymptotic  range.  Based  upon  this  observation,  it  is  entirely  possible  that  the 
grid  G2  does  lie  within  the  asymptotic  range  indicating  formal  convergence  of  the  tidal  solution. 
Note  that  the  addition  of  finite  amplitude  nonlinearities  results  in  a  periodic  semidiurnal  spike  in 
the  overprediction  error. 


7.2  Sensitivity  of  Nonlinearities  to  Mesh  Resolution 

Richardson  error  estimates  are  also  used  to  examine  the  relation  between  mesh  resolution,  the 
inclusion  of  various  nonlinear  terms,  and  the  model-computed  tidal  solution.  In  comparing  tidal 
elevation  predictions  derived  from  simulations  that  include  or  exclude  the  finite  amplitude  terms, 
Fig.  32a  and  b  respectively,  finite  amplitude  terms  appear  to  increase  the  overprediction  of  eleva¬ 
tion  as  noted  previously,  suggesting  that  finer  resolution  is  required  to  adequately  capture  the 
added  complexity  and  nonlinearity  associated  with  these  terms.  This  is  not  a  surprising  result,  since 
the  finite  amplitude  contribution  is  based  on  sea  surface  gradients  that  inherently  require  two  points 
for  computation.  The  increased  overprediction  may  also  be  due  to  a  greater  sensitivity  of  the 
nonlinear  dynamics  to  bathymetry.  Resolving  changes  in  the  bathymetry  gradient  requires  higher 
nodal  densities,  increasing  mesh  resolution. 
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FINE  GRID  COMPARISON:  G4  vs.  G9 
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EXTREME  ERROR  (cm) 
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Fig.  31 — Grid  comparison  simulations,  (a)  nonlinear  bottom  only  and  (b)  nonlinear  bottom  friction  and 

nonlinear  finite  amplitude  terms 


In  a  further  comparison,  Richardson-based  error  estimates  for  grid  G9  for  tidal  solutions  that 
have  either  the  convective  acceleration  terms  (space  derivatives)  active  or  inactive  are  displayed 
in  Fig.  32a  and  b.  Inclusion  of  the  convective  terms  not  only  leads  to  additional  peaks  in  the 
overprediction  error  curve  in  time,  but  the  underprediction  of  elevation  is  also  enhanced.  Again, 
finer  grid  resolution  is  necessary  to  adequately  capture  contributions  from  nonlinear  advective 
interactions  of  the  tide,  though  the  specific  cause  of  the  underprediction  is  uncertain.  Using  the 
finest  grid,  G2,  the  inclusion  of  convective  acceleration  terms  (space  derivatives)  leads  to  model 
instability.  One  shortcoming  of  the  ADCIRC-2DDI  model  in  its  current  state  is  the  spatial  discretization 
of  the  convective  terms,  coupled  with  an  explicit  treatment  in  time  of  nonlinearities  associated  with 
these  terms.  This  aspect  of  the  model  is  currently  under  investigation  with  improvements  expected 
in  subsequent  versions  of  the  code  (greater  than  v31_06).  Note  that  normalized  errors  for  all  cases 
are  relatively  small  and  uniform  throughout  the  simulation  period.  Though  errors  are  not  completely 
eliminated,  the  goal  of  minimizing  errors  associated  with  the  mesh  resolution  is  achieved. 
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EXTREME  ERROR  (cm) 
NORMALIZED  ERROR  (%) 


Fig.  32  —  Richardson-based  error  simulations,  (a)  fully  nonlinear,  (b)  nonlinear  bottom  friction/nonlinear 
convective  acceleration  terms,  and  (c)  nonlinear  bottom  friction  only 
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The  resolution  requirements  for  tidal  elevations  in  the  North  Sea/English  Channel  as  determined 
from  this  grid  convergence  study  indicate  that  nonlinear  dynamics  in  the  shallow  water  is  important 
and  requires  significantly  more  spatial  resolution.  Tidal  computations  are  highly  sensitive  to  bathymetric 
changes  and  coastline  geometry  that  occur  at  small  spatial  scales.  For  the  meshes  implemented, 
even  though  grid  G4  (whose  resolution  ranges  from  4-8  km)  is  not  within  the  asymptotic  range, 
it  does  seem  to  adequately  capture  smaller  scale  tidal  phenomena.  Concurrently,  the  mesh  does  not 
require  excessive  computer  resources.  For  any  study  of  tidal  dynamics,  it  is  extremely  important  to 
conduct  a  grid  convergence  check  to  minimize  influences  of  the  spatial  discretization  on  the  computed 
solution. 


8.0  CONCLUSIONS 

There  is  a  movement  in  recent  times  toward  the  application  of  coastal  ocean  models  that 
employ  FE  modeling  techniques.  Significant  advantages  are  associated  with  the  FE  discrete  form 
of  the  governing  equations.  The  most  obvious  of  these  advantages  is  the  tremendous  grid  flexibility 
afforded  by  FEs.  Variable  mesh  resolution  permits  mesh  refinement  in  specified  regions  to  resolve 
sharp  gradients  in  flow  or  bathymetric  features  and  to  capture  the  tortuous  detail  of  the  shoreline. 
At  the  same  time,  coarse  mesh  spacing  remains  in  deeper  waters  where  changes  in  the  ocean 
dynamics  are  known  to  occur  more  slowly  or  on  larger  scales. 

An  in-depth  examination  of  the  ADCIRC-2DDI  FE,  hydrodynamic  model  and  its  application  to 
tidal  prediction  in  coastal  regions  is  the  focus  of  this  report.  An  overview  of  the  basic  premise  of 
FE  approximations  and  the  mathematical  development  of  the  FE  equations  implemented  within  the 
ADCIRC  model  are  presented.  Utilization  of  the  ADCIRC-2DDI  model  in  no  way  requires  a  complete 
understanding  of  the  mathematical  foundation  of  the  FE  equation  development.  Rather,  the  user  is 
required  to  understand  the  relationship  between  the  model  assumptions,  mesh  resolution,  domain, 
boundary  forcing,  and  the  model-computed  fields  to  adequately  assess  the  computed  tidal  predictions. 

Using  the  TFF  benchmark  case,  the  ADCIRC-2DDI  model  is  validated  in  the  simulation  of  both 
tidal  elevations  and  currents.  Agreement  of  the  “base  simulation”  elevation  time  series  with 
observations  is  excellent  at  10  of  the  11  coastal  elevation  gauge  locations.  ADCIRC-2DDI  has 
demonstrated  accurate  prediction  of  the  tidal  elevations  at  these  11  coastal  stations  (during  a 
particular  25-h  time  period)  that  is  superior  to  several  other  tidal  models.  However,  due  to  the 
limitations  of  model-model  comparisons  using  previously  published  results,  a  conclusion  that  ADCIRC 
is  a  better  model  is  not  warranted.  As  a  side  note,  a  true  model-data  comparison  only  occurs  when 
the  modeled  and  observed  time  series  are  composed  of  identical  tidal  frequencies.  Considerable 
discussion  in  the  text  is  devoted  to  this  concept  within  the  context  of  the  North  Sea/English 
Channel  validation. 

Considering  the  quality  of  data,  the  agreement  of  the  velocity  hindcasts  is  also  quite  good, 
particularly  with  respect  to  current  directions.  For  speed  comparisons,  agreement  is  best  at 
velocity  gauges  placed  in  a  vertical  location  where  actual  speeds  are  expected  to  be  similar  to 
depth-averaged  speeds. 

Through  harmonic  analysis  and  recomposition  of  model  results,  the  11  tidal  constituents  included 
in  the  TFF  data  are  found  to  be  inadequate  for  a  complete  representation  of  the  actual  tidal  signal. 
From  model  solutions,  it  is  clear  that  several  additional  nonlinear  tides  make  significant  contributions 
to  the  elevation  in  near-coastal  areas  where  the  gauges  are  located. 
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The  sensitivity  analyses  undertaken  here  illuminates  the  influence  of  nonlinearities,  forcing, 
and  several  model  parameters  on  the  model-computed  response.  Modular  construction  of  the 
ADCIRC-2DDI  model  makes  such  analyses  possible  without  requiring  user  modification  of  the  code. 
Several  conclusions  and  recommendations  based  on  these  analyses  are  made  regarding  the  use  of 
ADCIRC-2DDI.  First,  due  to  its  physical  basis  and  mass  conservation  properties,  the  nonlinear 
parameterization  of  bottom  stress  offers  measurable  advantages  over  use  of  a  linear  bottom  friction 
formulation.  Secondly,  the  inclusion  of  other  nonlinear  terms  in  the  GWCE  (such  as  the  finite 
amplitude  or  convective  terms)  results  in  noticeable,  though  not  dominant,  differences  in  the  computed 
model  solution.  Bottom  friction  remains  the  single  most  important  mechanism  affecting  tidal  elevations 
and  circulation  in  the  coastal  environment.  A  fully  nonlinear  simulation  in  which  all  nonlinear 
terms  contribute  to  the  model  dynamics  offers  a  more  complete  representation  of  tidal  physics, 
especially  with  respect  to  internal  nonlinear  tide  generation. 

For  the  North  Sea/English  Channel  simulations,  a  fully  nonlinear  model  response  produces 
more  accurate  model-to-data  comparisons  with  respect  to  tidal  elevation.  Care  must  be  taken  when 
exercising  the  model  in  a  fully  nonlinear  mode  to  maintain  an  appropriately  fine  mesh  resolution 
and  avoid  the  drying  of  elements  through  a  minimum  depth  specification  if  wetting/drying  is 
inactive. 

This  sensitivity  study  confirms  that  nonlinear  tides  on  the  shelf  can  be  important  and  should 
be  included  in  tidal  boundary  forcing  when  open  boundaries  are  located  on  the  continental  shelf. 
Alternatively,  utilizing  a  larger  domain  locates  the  open  boundaries  in  deeper  water  where  the 
nonlinearities  are  of  little  consequence.  For  domains  larger  than  the  North  Sea/English  Channel 
region  studied,  O  (1000  km  x  1000  km),  tidal  potential  forcing  and  a  variable  Coriolis  forcing  will 
likely  have  a  greater  significance.  The  range  of  latitudes  spanned  by  the  domain  strongly  influence 
the  importance  of  these  mechanisms.  In  the  context  of  the  southern  North  Sea/English  Channel 
case,  inclusion  of  horizontal  mixing  in  the  form  of  lateral  eddy  viscosity  is  of  little  consequence. 
If  the  environment  is  advection  dominated  or  small-scale  nearshore  flows  are  of  interest,  the  horizontal 
mixing  is  expected  to  increase  in  importance. 

Finally,  a  brief  examination  of  mesh  resolution  effects  on  the  computed  tidal  solution  are 
undertaken  in  the  context  of  a  grid  convergence  analysis.  Though  grids  having  9  km  (G9)  and  4  km 
(G4)  as  minimum  nodal  spacings  do  not  fall  in  the  asymptotic  range,  small  differences  between  the 
Richardson-based  error  estimate  for  grid  G4  suggests  that  a  grid  having  a  minimum  resolution  of 
2  km  (G2)  is  at  or  very  near  the  asymptotic  range.  Most  of  the  significant  errors  resulting  from 
insufficient  grid  resolution  are  thought  to  be  due  to  the  misrepresentation  of  large  velocity  gradients 
induced  by  the  irregular  coastline.  Model  nonlinearities  further  magnify  such  errors. 

The  primary  purpose  for  this  study  has  been  to  assess  the  capability  of  the  ADCIRC-2DDI 
barotropic,  2D,  FE  hydrodynamic  model  to  predict  tidal  dynamics  in  coastal  waters.  A  comprehensive 
series  of  experiments  testing  sensitivity  of  the  model  have  provided  considerable  insight,  not  only 
to  the  application  of  the  model  itself,  but  also  to  the  behavior  of  tides  in  a  coastal  environment.  The 
hope  is  that  experiences  gained  from  this  work  will  be  carried  over  into  future  tidal  applications 
leading  to  an  improved  predictive  capability. 
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Appendix  A 

RAW  MODEL-DATA  ELEVATION  COMPARISONS  FOR  THE  BASE  SIMULATION 


o  OBSERVATIONS 
- BASE  SIMULATION 


Fig.  A1  —  Base  simulation  tidal  elevations,  tidal  gauge  at  St.  Malo 
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Fig.  A2  —  Base  simulation  tidal  elevations,  tidal  gauge  at  Cherbourg 
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o  OBSERVATIONS 
-  BASE  SIMULATION 


Base  simulation  tidal  elevations,  tidal 
gauge  at  Dieppe 


Fig.  A4  —  Base  simulation  tidal  elevations,  tidal 
gauge  at  Boulogne 
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o  OBSERVATIONS 
- BASE  SIMULATION 


Fig.  A9  —  Base  simulation  tidal  elevations,  tidal 
gauge  at  Dover 
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Fig.  A10  —  Base  simulation  tidal  elevations,  tidal 
gauge  at  Christchurch 
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Fig.  All  — Base  simulation  tidal  elevations,  tidal 
gauge  at  Lowestoft 
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Appendix  B 


“FILTERED”  MODEL-DATA  ELEVATION  COMPARISONS 
FOR  THE  BASE  SIMULATION 


o  OBSERVED 

- RAW  MODEL  OUTPUT 

RECOMPOSED 

Fig.  B1  —  Base  simulation  tidal  elevations,  11  TFF 
constituents  only,  tidal  gauge  at  St.  Malo 


Fig.  B2  —  Base  simulation  tidal  elevations,  11  TFF 
constituents  only,  tidal  gauge  at  Cherbourg 
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DAY  (March  1976) 
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o  OBSERVED 

- RAW  MODEL  OUTPUT 

RECOMPOSED 


Fig.  B4  —  Base  simulation  tidal  elevations,  11  TFF 
constituents  only,  tidal  gauge  at  Boulogne 
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Fig.  B5  —  Base  simulation  tidal  elevations,  11  TFF 
constituents  only,  tidal  gauge  at  Calais 
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o  OBSERVED 
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Fig.  B9  —  Base  simulation  tidal  elevations,  11  TFF 
constituents  only,  tidal  gauge  at  Dover 
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Fig.  BIO  —  Base  simulation  tidal  elevations,  12 
TFF  constituents  only,  tidal  gauge  at  Christchurch 
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Fig.  Bll  —  Base  simulation  tidal  elevations,  11 
TFF  constituents  only,  tidal  gauge  at  Lowestoft 
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Fig.  C3  —  Base  simulation  constituent  amplitudes, 
tidal  gauge  at  Dieppe 
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Fig.  C4  —  Base  simulation  constituent  amplitudes, 
tidal  gauge  at  Boulogne 
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Fig.  C5  —  Base  simulation  constituent  amplitudes, 
tidal  gauge  at  Calais 
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Fig.  C12  —  Base  simulation  constituent  phases, 
tidal  gauge  at  St.  Malo 
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Fig.  C13  —  Base  simulation  constituent  phases, 
tidal  gauge  at  Cherbourg 
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Fig.  C14  —  Base  simulation  constituent  phases, 
tidal  gauge  at  Dieppe 
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Fig.  C15  —  Base  simulation  constituent  phases, 
tidal  gauge  at  Boulogne 
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Fig.  C16  —  Base  simulation  constituent  phases, 
tidal  gauge  at  Calais 
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Fig.  C17  —  Base  simulation  constituent  phases, 
tidal  gauge  at  Zeebrugge 
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Fig.  C24  —  Base  simulation  constituent  errors  in  the 
complex  plane,  tidal  gauge  at  Cherbourg 
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Fig.  C25  —  Base  simulation  constituent  errors  in  the 
complex  plane,  tidal  gauge  at  Dieppe 
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Fig.  C26  —  Base  simulation  constituent  errors  in  the 
complex  plane,  tidal  gauge  at  Boulogne 
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Fig.  C32  —  Base  simulation  constituent  errors  in 
the  complex  plane,  tidal  gauge  at  Christchurch 
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Fig.  C33  —  Base  simulation  constituent  errors  in 
the  complex  plane,  tidal  gauge  at  Lowestoft 
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Fig.  D3  —  Base  simulation  speeds,  velocity 
gauge  3 
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Fig.  D4  —  Base  simulation  speeds,  velocity 
gauge  4 
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Fig.  D5  —  Base  simulation  speeds,  velocity 
gauge  5 
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Fig.  D6  —  Base  simulation  speeds,  velocity 
gauge  6 
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Fig.  D7  —  Base  simulation  speeds,  velocity 
gauge  7 
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Fig.  D8  —  Base  simulation  speeds,  velocity 
gauge  8 
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